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A variational approach for constructing an effective particle description of the low-energy physics 
of one-dimensional quantum spin chains is presented. Based on the matrix product state formalism, 
we compute the one- and two-particle excitations as eigenstates of the full microscopic Hamiltonian. 
We interpret the excitations as particles on a strongly-correlated background with non-trivial dis¬ 
persion relations, spectral weights and two-particle S matrices. Based on this information, we show 
how to describe a finite density of excitations as an interacting gas of bosons, using its approxi¬ 
mate integrability at low densities. We apply our framework to the Heisenberg antiferromagnetic 
ladder: we compute the elementary excitation spectrum and the magnon-magnon S matrix, study 
the formation of bound states and determine both static and dynamic properties of the magnetized 
ladder. 
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I. INTRODUCTION 

Finding the ground state of strongly-correlated quan¬ 
tum many-body systems poses one of the main challenges 
of contemporary condensed matter physics. The physical 
properties of these systems, however, are not determined 
by the ground state but rather by the low-lying, elemen¬ 
tary excitations relative to this ground state. In contrast 
to the strongly-correlated ground state, these elementary 
excitations are of a particularly simple character: in most 
cases they can be treated as a collection of independent, 
weakly-interacting particles living on non-trivial vacuum 
stat^^l. 

In condensed matter theory these “quasi-particles” 
are typically defined starting from some non-interacting 
limit. In Fermi liquid theory for interacting electron 
system^EI - the most prominent example of this ap¬ 
proach - the quasi-particles are defined in the free- 
electron system, but remain well-defined modes when 
turning on the interactions. The effect of a finite life¬ 
time and quasi-particle interactions can be treated in per¬ 
turbation theory. In strongly-correlated lattice systems, 
however, there is typically no obvious way to start from 
a non-interacting theory to define the quasi-particles 
that determine the system’s properties (notable counter¬ 
examples in one dimension include integrable systems 5 
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and continuous unitary transformations®). The varia¬ 
tional approach, which we will advocate in this paper, is 
orthogonal to the perturbative approach by starting from 
the strongly-correlated ground state and finding the low- 
lying excitations of the interacting system variationally. 
As exact eigenstates, these excitations have an infinite 
lifetime, but a priori it is not clear that they should have 
a local, particle-like nature. 

In relativistic quantum held theory, a picture of lo¬ 
calized elementary excitations on top of a strongly- 
correlated vacuum has been formulated in a rigorous 
fashion. Haag-Ruelle scattering theorjf? does indeed con¬ 
struct a many-particle Fock space by acting with local 
operators on the vacuum and even defines an S matrix 
describing the interactions between these particles. Be¬ 
cause this formalism depends heavily on Lorentz invari¬ 
ance, there is a priori no straightforward translation to 
lattice systems. Indeed, on the lattice there are fewer 
restrictions on the spectrum: different elementary ex¬ 
citation branches typically have different characteristic 
velocities and are not bound to be stable in the whole 
Brillouin zone - a typical spectrum is shown in Fig. [I] 

Recently though, it was realized that by using Lieb- 
Robinson bounds® as the soft lattice analog of strict 
causality in relativistic QFT, the locality of elementary 
excitations can be established in a rigorous way. More 
specifically, it was shown in Ref. IS] that an elementary 
excitation that lives on an isolated branch in the energy- 
momentum spectrum and has a finite overlap with an ar¬ 
bitrary local operator, can be created out of the ground 
state by the action of a momentum superposition of a 
local operator (to an exponential precision in the size of 
the support of this operator). In Ref. jTO] the scattering 
problem of these particle excitations was formulated by 
translating the Haag-Ruelle formalism to the lattice set¬ 
ting. 

These theoretical developments provide a motivation 
for the variational approach towards a particle picture 
of the low-energy physics of lattice systems. Indeed, by 
making use of the fact that gapped excitations should 
be local, it might prove possible to describe them with 
only a small number of variational parameters. This is 
the idea of the single-mode approximation, or Feynman- 
Bijl ansatz, pioneered by Feynman in his study on liq¬ 
uid helium 11 21 and later successfully applied to quan¬ 
tum spin system^SHUl. 

Although providing qualitative 
insight into the nature of the elementary excitations, the 
single-mode approximation is often too crude as a vari¬ 
ational ansatz to obtain quantitative results on the low- 
lying spectrum of generic spin chains. 

Indeed, constructing a variational ansatz for excita¬ 
tions with quantitative accuracy requires both an accu¬ 
rate parametrization of the ground state and a system¬ 
atic way to change this ground state locally. For one¬ 
dimensional systems, the framework of matrix product 
states 16 17 (MPS) has proven to meet both requirements. 
The ground state of one-dimensional quantum spin sys¬ 
tems can indeed be efficiently parametrized by the class of 
MpcjH the success of the density matrix renormaliza¬ 
tion groupPHis based on MPS serving as the class of states 
over which it optimize d — I— ] The defining characteristic 


of MPS - or tensor network states in general - is the pres¬ 
ence of a “virtual” level that carries the (quantum) corre¬ 
lations in the many-body wave function. By acting both 
on the physical and the virtual level, a variational ansatz 
for elementary excitations on an MPS ground state was 
introduced in Refs. [23] and |2U The ansatz was used for 
calculating dispersion relations and dy namical correla¬ 
tion fu nction s of quantum spin chain d 25 '' 27 l, quantum field 
theoried 2 ® 221 and gauge theoried 3 ®. 

A more general understanding of these efforts is ob¬ 
tained by realizing that the low-energy dynamics corre¬ 
spond to small variations around the variational ground 
state and are therefore not necessarily contained within 
the variational class itself. Indeed, for the smooth 
manifold of MPfpS it is the tangent space constructed 
around the MPS ground state that provides a natural 
parametrization of the low-energy dynamics. For ex¬ 
ample, the best approximation to time evolution within 
the MPS manifold can be obtained by projecting the 
Schrodinger equation into the MPS tangent bundle, 
according to the time-dependent variational principle 
(TDVP) 32 33 : Similarly, the ansatz for an elementary 
excitation corresponds exactly to a vector in the tangent 
space around the MPS ground state. These ideas can be 
grouped under the concept of post-MPS methods 2 ® as an 
alternative for the standard MPS algorithms for tackling 
the low-energy dynamics around an MPS ground state. 

The crucial next step in this approach - after the con¬ 
struction of single-particle excitations - consists of study¬ 
ing the interactions between these particles and, more 
specifically, computing the two-particle S matrix 3 ®. This 
information can then be used as the input for the “ap¬ 
proximate Bethe ansatz”!^ 3 ® (i.e. neglecting all three- 
particle scattering processes) in order to provide a first- 
quantized description of a finite density of excitations on 
top of the strongly-correlated vacuum. 

These developments should eventually lead to the 
ab initio construction of an effective second-quantized 
Hamiltonian, acting in a Fock space of interacting par¬ 
ticles. In contrast to standard effective field theory con¬ 
structions, the variational approach would automatically 
incorporate all symmetries and correlations of the vac¬ 
uum state on which these particles live without relying 
on phenomenological considerations. 

In this paper we further elaborate on the framework 
that we introduced in Ref. ]3H In Sec. El we show how to 
construct one- and two-particle excitations on an MPS 
vacuum state. We formulate a definition of the S matrix 
based on the form of the two-particle wave function and 
prove that it is unitary. Finally, we construct the pro¬ 
jector on the one- and two-particle subspace which shows 
up in the spectral representation of dynamical correlation 
functions. In Sec. m we take a step back and show that 
the S matrix as defined in Sec. |TT| corresponds to the one 
that shows up in standard dynamical scattering theory. 
Next we elaborate on the approximate Bethe ansatz as 
a way of dealing with a finite density of excitations in 
a first-quantized many-particle formalism. In Sec. |IV| we 
apply our variational method to study the Heisenberg an¬ 
tiferromagnetic two-leg ladder. We compute the elemen¬ 
tary excitation spectrum, the two-particle S matrix and 
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FIG. 1. A typical momentum-energy excitation spectrum of 
a one-dimensional lattice system. We have depicted three el¬ 
ementary (one-particle) excitations (full lines) and the many- 
particle continuum (grey). Both «2 and 03 are stable in the 
whole Brillouin zone; the latter remains stable even inside the 
continuum, possibly because it cannot decay in a two-particle 
state through symmetry constraints. Particle ai becomes un¬ 
stable upon entering the continuum, so that it ceases to be 
a one-particle excitation (cannot be created by a local opera¬ 
tor). 


one- and two-particle contributions to dynamical corre¬ 
lation functions. Afterwards, we apply the approximate 
Bethe ansatz to the magnetization process, at zero and 
finite temperature, and compute both thermodynamic 
properties and correlation functions of the magnetized 
ladder. In the last section, we provide an overview of 
some interesting extensions of our framework and give an 
outlook towards the construction of effective held theo¬ 
ries in second quantization. 


II. CONSTRUCTING SCATTERING STATES 


In this section we construct variational one- and two- 
particle states on an MPS background®. Based on the 
form of these wave functions, we define the S matrix and 
introduce the projectors on the one- and two-particle sub¬ 
spaces (i.e. the low-energy subspace) which can be used 
to compute the low-energy part of dynamical correlation 
functions. Note that, while the complete framework is 
presented in the main body, technical details and long 
equations are taken up in the appendix. A short remark 
on notation is also in order. Vectors of any length will 
be denoted in boldface, whereas matrices will use a sans 
serif font. Vector entries will be referred to using a su¬ 
perscript (in which case the boldface will be dropped), 
whereas subscripts of a boldface vector typically refer to 
a label of a set of vectors, such as a basis. The only ex¬ 
ception to these rules is that physical states are denoted 
using Dirac’s bra-ket notation and the matrices appear¬ 
ing in the definition of the matrix product state (which 
can also be interpreted as rank three tensors) are typeset 
using the normal serif type (italic). 


i 
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FIG. 2. Graphical representation of an MPS ground state. 
The circles represent the (D x d x D)-dimensional tensor A: 
every outgoing leg corresponds to a tensor index. Whenever 
two legs are connected, this corresponds to a contraction of 
the two indices. In the MPS all virtual indices are contracted, 
while the physical indices correspond to the physical degrees 
of freedom in the MPS wave function |l]). The matrix product 
structure contains the (quantum) correlations of this ground 
state. 


A. Ground state 


Consider a one-dimensional quantum spin system with 
local dimension d in the thermodynamic limit, described 
by a local and translation invariant Hamiltonian. While 
in no way crucial, we restrict to nearest-neighbour Hamil¬ 
tonians, i.e. H = ^, !gZ h n ,n+i) for reasons of simplic¬ 
ity. We furthermore assume that the translation invari¬ 
ant ground state of this system (we restrict to a unique 
ground state, see Sec.[V]for extensions) can be accurately 
described by an injective uniform matrix product state 
(uMPs^mni 


\m) = E v l 

M=i 


11 1 

m 


v r IM>» 


(i) 


where A s is a set of D x D matrices for s = 1,..., d, or, 
equivalently, A can be interpreted as a D x d x D tensor; 
v^ L and vr are D-dimensional boundary vectors. In the 
thermodynamic limit, all physical observables are inde¬ 
pendent of these boundary vector^®, so that the tensor 
A provides a complete description of the ground state 

The set of injective MPS of a certain bond dimension 
constitute a complex manifold®. Finding the best ap¬ 
proximation of the ground state of a certain Hamilto¬ 
nian within this manifold can be achieved using different 
algorithm^® - in_our simulations we will always use the 
TDVP algorithnPmi 


B. One-particle excitations 

This ground state serves as our vacuum, on top of 
which we will build localized, particle-like excitations. A 
first guess for the wave function of an elementary excita¬ 
tion with momentum k is the single-mode approximation 

|Ssma(k)> = 5> im 6 n |tf[.4]>, (2) 

n 

where O n is an operator acting at site n. The c hoice 
of operator can be inspired by physical intuitioii 1 d l 15 l 41 l 
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FIG. 3. Graphical representation of the one-particle excita¬ 
tion ansatz. The ground state tensor A is changed at site n 
into a new tensor B (square) and a momentum superposition 
is taken. The matrix product structure allows that the tensor 
B can change the ground state over a finite distance. 


or determined by numerical optimization 42 . Though pro¬ 
viding some qualitative insight into elementary excitation 
spectra, this ansatz is typically not a good quantitative 
approximation for the true wave function of the excita¬ 
tion. Systematically improving on this would ask for the 
introduction of bigger local operators O n . It was indeed 
proven^ that, in the case of an isolated excitation branch, 
the exact wave function can be arbitrary well approxi¬ 
mated in this way. More specifically, it was shown that 
the localized nature of an excitation depends on the gap 
to the nearest eigenvalue of the Hamiltonian in the same 
momentum sector. 

Within the framework of matrix product states, it is 
possible to construct a variational ansatz that is able to 
capture the localized nature of the excitation by directly 
modifying the local tensors. Indeed, instead of only op¬ 
erating on the physical level, we can change one ground 
state tensor A s wi th a n ew tensor B s and take a momen¬ 
tum superpositiorpil23I25| 


|$ K [ J B])=^e i -^4 

n { s } 


11 '• 

_m<n 


X 


II V 

_m>n 


wilM)- 


(3) 


Through the virtual level of the MPS, this ansatz is able 
to perturb the ground state over a finite length deter¬ 
mined by the bond dimension D. All variational freedom 
of this ansatz is contained within the tensor B s . As the 
parametrization of the state ^ is linear in the elements 
of B s , variationally optimizing amounts to solving the 
Rayleigh-Ritz problem 


Hcff,i P (re)M = Am (4) 

with H e fy lp («) the momentum dependent effective one- 
particle Hamiltonian and vector u containing the coeffi¬ 
cients u l to expand tensor B in the state ([3]) with respect 
to a suitably chosen basis {B^,i = 1,..., (d — 1 )D 2 }. 
We refer to Appendix |A 2| for details on how to calculate 

Heff,lp ■ 

Upon solving the eigenvalue problem in Eq. Q , we ob¬ 
tain a set of ( d — 1 )D 2 eigenvalues for every momentum k. 
Some of those eigenvalues A a (/%) offer a good approxima¬ 
tion to the exact dispersion relations of the elementary 
excitations, i.e. the isolated branches in the spectrum 
of the Hamiltonian. Moreover, from the corresponding 
eigenvectors m q (k) we obtain an explicit expression for 


the wave function of the elementary excitations by insert¬ 
ing the tensors B a (n) = Y^i u a( K )B(i) hr Eq. This 
expression can be used to calculate the spectral weights 
of the excitations and, consequently, the one-particle con¬ 
tribution to dynamical correlation functions. 

Other eigenvalues obtained from Eq. Q will fall in the 
continuous part of the spectrum of the Hamiltonian, i.e. 
in the set of scattering states. Scattering states cannot 
be described by a single local perturbation, so we expect 
the ansatz ([3]) to fail. In fact, instead of a scattering 
state, the variational optimization will create a localized 
wave packet of two-particle states within some energy 
range. Obviously, the variational eigenstates of the form 
in Eq. (J3| will not provide a good approximation to the 
exact scattering eigenstates of the full Hamiltonian. A 
more appropriate variational ansatz for two-particle scat¬ 
tering states is discussed in the remainder of this section. 

Remark that we have so far not discussed the case of 
bound states. When defining (quasi-) particles along a 
path of Hamiltonians using e.g. perturbation theory or 
continuous unitary transformation^, bound states can 
be identified with isolated eigenstates emerging from a 
multi-particle continuum along the path. In our vari¬ 
ational framework, we consider one particular Hamilto¬ 
nian which is not necessarily related to a one-parameter 
family. All isolated branches in the spectrum are equally 
elementary (see Ref. l43l for the analogous result in QFT). 
While there might be quantum numbers that indicate 
the “history” of an elementary excitation along a path of 
Hamiltonians, there is typically no particle number sym¬ 
metry to make the interpretation of bound states unam¬ 
biguous. On a related note, elementary excitations are by 
this definition exact eigenstates of the Hamiltonian and 
therefore have an infinite life time. We cannot and do 
not target resonances within the continuous part of the 
spectrum. Therefore, the Hamiltonian does not contain 
interactions that link the one-particle sector with higher 
particle states. 

As mentioned previously, the spectrum of general 
quantum spin chains can be very complex. Within cer¬ 
tain regions of the Brillouin zone, the energy of elemen¬ 
tary excitations can fall within the continuum (this typ¬ 
ically requires a quantum number that protects them 
against decay) or there might be no elementary excita¬ 
tions at all (e.g. around momentum zero in the spin-1 
Heisenberg antiferromagnet). We therefore need a way 
to determine which variational eigenvalues of Eq. Q cor¬ 
respond to elementary excitations and therefore offer a 
good approximation to actual eigenstates of the Hamilto¬ 
nian. Upon enlarging the variational one-particle space, 
e.g. by increasing the bond dimension or the spatial 
support of the local perturbation, eigenvalues that cor¬ 
respond to elementary excitations will converge quickly 
(related to the gap to the nearest eigenvalue) and re¬ 
main at a fixed position. Eigenvalues in the continu¬ 
ous part of the exact spectrum, on the other hand, will 
not really converge and several new eigenvalues will ap¬ 
pear in those regions. A more quantitative way to as¬ 
sess how well an exact eigenstate is approximated con¬ 
sists of calculating the variance of the HamiltoniarP^, i.e. 
(<fr K [i?]| (H — A(k)) 2 |4> k [H]). For elementary excitations, 
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these variances should be small (see Sec. IV A for numer¬ 
ical values). For the other solutions of the one-particle 
problem Q, which correspond to scattering states, the 
variance should be larger. For a typical gapped system, 
the difference will be some orders of magnitude. Con¬ 
sequently, this quantity allows for the identification of 
one-particle states, even within higher-particle bands and 
without exploiting symmetries. 

Note finally that, without Galilean invariance on the 
lattice, the tensor B a (n), which describes the parti¬ 
cle a on a dispersion branch A„(k), is momentum de¬ 
pendent. On the other hand, we expect that for a 
well-defined particle in a certain momentum range this 
momentum dependence is not too strong. Indeed, it 
turns out that by a suitable choice of the basis tensors 
= 1,..., (d — 1)-D 2 }, we can fully capture B a (n) 
for all elementary excitations a and for all momenta k 
in the span of just a small number t -C (d — 1 )D 2 basis 
vectors {B^, i = 1,..., t}. Although more sophisticated 
optimization strategies should be possible, we construct 
this reduced basis from a number of B' s at different mo¬ 
menta. This reduced basis will be important for solving 
the scattering problem in the next sections. 


C. Variational ansatz for two-particle states 

In the previous section it became clear that we need 
another ansatz to capture the delocalized nature of a two- 


particle state. We will start from a one-particle spectrum 
consisting of a number of different types of particles, la¬ 
belled by a, with dispersion relations A q (k). In the ther¬ 
modynamic limit, constructing the two-particle spectrum 
is trivial: the momentum and energy are the sum of the 
individual momenta and energies of the two particles 2 . 
The two-particle wave function, however, depends on the 
particle interaction. The interactions, which depend on 
both the Hamiltonian and the ground state correlations, 
are reflected in the wave function in two ways: (i) the 
asymptotic wave function has different terms, with the S 
matrix elements as the relative coefficients, and (ii) the 
local part of the wave function. 


In order to capture both we introduce the following 
ansatz for describing states with two localized, particle¬ 
like excitations with total momentum K 


+oo L n 

i T w> = EE^)ix*>)> (5) 

n=0 j—1 


where the basis states are 


B 


+oo d 

\ X K,j(n = 0 )) = E * lKni E II ^ 

m =—oo {s}=l [_ra<ni 

+oo d 

\XK, Uuh) (n>0))= E * lKni E v l U AS 

m=—oo {s}=l [_m<ni 


U) 


We collect the variational coefficients either in one half¬ 
infinite vector C with C J,rl = c^(n) or using the finite 
vectors c(n) with entries {c J (n), j = 1,..., L n } for every 

n = 0,1,_ Here, we have L 0 = (d— 1 )D 2 and L n>0 = 

[(d — 1)Z? 2 ] 2 . Note that the sum in Eq. © only runs over 
values n > 0, because a sum over all integers would result 
in an overcomplete basis. 

Already at this point, we will reduce the number of 
variational parameters to keep the problem tractable. 
The terms with n = 0 (corresponding to the basis vec¬ 
tors in Eq. ©) are designed to capture the situation 
where the two particles are close together. No infor¬ 
mation on how this part should look like is a priori 
available, so we keep all variational parameters c° (0), 
j = 1 ,... ,L 0 = D 2 (d — 1). The terms with n > 0 cor¬ 
responding to the basis vectors in Eq. 0 represent the 
situation where the particles are separated. We know 
that, as n —> oo, the particles decouple and we should 
obtain a combination of one-particle solutions. With this 
in mind, we restrict the range of j\ and jz to the first £ 


II A ‘ 

_m>n i 




( 6 ) 


B r\ 

Ui) 

n • 4 - 

(J2) 

n - 1 *" 


_ni<ra<ni+n 


_ra>ni+n 


VR |{ s }) • (?) 


basis tensors {B^,i = 1,...,£}, which were chosen so 
as to capture the momentum dependent solutions of the 
one-particle problem. Consequently, the number of basis 
states of Eq. 0 for n > 0 satisfies L n = £ 2 , which we 
will henceforth denote as just L. 

This might seem like a big approximation for n small: 
when the two particles approach the wave functions 
might begin to deform, so that the B tensors of the one- 
particle problem no longer apply. Note, however, that 
the local (n = 0) and non-local (n > 0) part are not or¬ 
thogonal, so that the local part is able to correct for the 
part of the non-local wave function where the one-particle 
description is no longer valid. 

As the state 0 is linear in its variational parameters 
C, optimizing the energy amounts to solving a general¬ 
ized eigenvalue problem 


H eff C = wN off C 


( 8 ) 
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FIG. 4. Graphical representation of the basis states 0 - The ground state is changed at two sites at a distance of n sites and 
a momentum superposition is taken (with the total momentum K). 


with to the total energy of the state and 

(H e s)n'j',nj = )l H \X],k[p)) (9) 

(NeflOra'jhnj = (Xj' : I<( n )| Xj,K(jl)) (10) 

two half-infinite matrices. They have a block matrix 
structure, where the submatrices are labelled by ( n',n ) 
and are of size L n i x L n . The computation of the matrix 
elements is quite involved and technical, so we refer to 
the appendix for the explicit formulas. 

Since the eigenvalue problem is still infinite, it can¬ 
not be diagonalized straightforwardly. Since we actually 
know the possible energies lj for a scattering state with 
total momentum K , we can also interpret Eq. § as an 
overdetermined system of linear equations for the coeffi¬ 
cients Ci ,n = c?(n). In the next two sections we will show 
how to reduce this problem to a finite linear equation. 


Hamiltonian H^ip, we already know a number of solu¬ 
tions to Eq. (13). Indeed, if we can find T combinations 
of two types of particles (a,/3) with individual momenta 
(/ci, K 2 ) such that K = K 1 +K 2 and lo = A a (/ci)+A | g(/C 2 ), 
then the polynomial eigenvalue problem will have 2 T so¬ 
lutions fi on the unit circle. These solutions take the form 
/i = e lK2 and the corresponding eigenvector is given by 


v = u a (n i)®up(K 2 ) (14) 

(in the case of degenerate eigenvalues we can take linear 
combinations of these eigenvectors that no longer have 
this product structure). Every combination is counted 
twice, because we can have particle a on the left and 
particle p on the right, and vice versa. 

Moreover, since Aj n = A_ m , the number of eigenvalues 
within and outside the unit circle should be equal. This 
allows for a classification of the eigenvalues /i as 


D. Asymptotic regime 

First we solve the problem in the asymptotic regime, 
where the two particles are completely decoupled. This 
regime corresponds to the limit n', n —» 00 , where the 
effective norm and Hamiltonian matrices, consisting of 
blocks of size L x L, take on a simple form. Indeed, if 
we properly normalize the basis states, the asymptotic 
form of the effective norm matrix reduces to the identity, 
while the effective Hamiltonian matrix is a repeating row 
of block matrices centred around the diagonal 

(Heff)n',n f A n _ n /, Tl : Tl / OO. (H) 

The blocks decrease exponentially as we go further from 
the diagonal, so we can, in order to solve the problem, 
consider them to be zero if \n— n'\ > M for some suitably 
chosen integer M. In this approximation, the coefficients 
c(n) obey 

M 

A m c( n + m)=ojc(n), n —> 00 . ( 12 ) 

m=—M 

We can reformulate this as a recurrence relation for the 
c(n) vectors and therefore look for elementary solutions 
of the form c(n) = \i n v. For fixed w, the solutions /i 
and v are now determined by the polynomial eigenvalue 
equation 

M 

y = OJV. (13) 

m=—M 

From the s pecial structure of the blocks A rn (see Ap- 
pendix|B~3|) and their relation to the one-particle effective 


I/aI 

< 1 

for 

i = 1,...,LM-T 


I/aI 

= 1 

for 

i = LM- r + 1 ,.. 

., LM + T 

I/aI 

> 1 

for 

i = LM + T + 1,.. 

., 2 LM. 


The last eigenvalues with modulus bigger than one are 
not physical (because the corresponding c(n) ~ 
yiels a non-normalizable state) and should be discarded. 
The 2T eigenvalues with modulus 1 are the oscillating 
modes discussed above; we will henceforth label them 
with 7 = 1 ,..., 2 T such that /r = e lK ' 1 ' (/t 7 being the 
momentum of the particle of the right) and the corre¬ 
sponding eigenvector is given by 

Vj = u ai (K - Kj) <g> u ^ (k 7 ). 

Finally, the first eigenvalues are exponentially decreas¬ 
ing and represent corrections when the excitations are 
close to each other. We henceforth denote them as e~ Xi 
with 5ft(Ai) > 0 for * = 1,..., LM — T and denote the 
corresponding eigenvectors as w t . 

With these solutions, we can represent the general 
asymptotic solution as 

LM-T 2T 

c(n) y p i e- Xin w i + yq 1 e iK ~< n v 1 . (15) 

i —1 7=1 

Of course, we still have to determine the coefficients 
{jAq 7 } by solving the local problem. 

E. Solving the full eigenvalue equation 

Since the energy uj was fixed when constructing the 
asymptotic solution, the generalized eigenvalue equation 
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is reduced to the linear equation 

(H e fi - wN cff )C = 0. 


We know that in the asymptotic regime this equation is 
fulfilled if and only if c(n) is of the form of Eq. (151. 
We will introduce the approximation that the elements 
for the effective Hamiltonian matrix [Eq. ([9])] and norm 
matrix [Eq. ( |10| ] have reached their asymptotic values 
when either n > M + N or n' > M + N, where N is 
a finite value and can be chosen sufficiently large. This 
implies that we can safely insert the asymptotic form for 
n > N in the wave function, which we can implement by 
rewriting the wave function as 


C = Zxx 


(16) 


where 


Z = 



{e Xin wi} 



The {e~ XiH Wi} and {e~ lK ~i n v 1 } are the vectors corre¬ 
sponding to the damped, resp. oscillating modes, while 
the identity matrix is inserted to leave open all parame¬ 
ters in c(n ) for n < N. The number of parameters in x is 
reduced to the finite value of D 2 {d — 1) + NL + LM + T. 

Since the equation is automatically fulfilled after M + 
N rows, we can reduce H c ff and N c ff to the first rows, so 
we end up with the following linear equation 


[H — wN] re d xZxi = 0 


(17) 


with 


[H — wN] red = 


/ 


0 

0 

... 0 

\ 



0 

0 

... 0 


(H 

- wN)ex 

A M 

0 

... 0 




Am-i 

Aa/ 

... 0 


V 


Ar 

a 2 

• ■ • Am 

/ 


We start from T linearly independent scattering eigen¬ 
states \Ti(K, lo)) (i = 1 ,...,T) at total momentum I\ 
and energy ui with asymptotic coefficients qi(K,u>). The 
asymptotic form of these eigenstates is thus a linear 
combination of all possible non-decaying solutions of the 
asymptotic problem: 


2 r 

| r i {K,u)) = Y l QUK,u) 

7=1 

(18) 

n>N j 


where the coefficients are obtained from solving the lo¬ 
cal problem. The number of eigenstates equals half the 
number of oscillating modes that appear in the linear 
combination. With every oscillating mode 7 we can as¬ 
sociate a function w 7 (k) giving the energy of this mode 
as a function of the momentum k 7 of the second particle 
at a fixed total momentum K. If 7 corresponds to the 
two-particle mode with particles a 7 and /3 7 , this function 
is given by w 7 (ft) = A a ^(K— k) + A^ 7 (k). The derivative 
of this function, which will prove of crucial importance, 
is w 7 (k) = A p (ft) — A ' a (K — ft). It can be interpreted 
as the difference in group velocity between the two parti¬ 
cles, i.e. the relative group velocity in the center of mass 
frame. 

Much like the proof of conservation of particle current 
in one-particle quantum mechanics, it can be shown that 
(see Appendix [C|, if (18) is to be the asymptotic form of 
an eigenstate, the coefficients q] (K. u>) should obey 


Ei^-)I 2 (^m)=o- (19) 


This equation can indeed be read as a form of conserva¬ 
tion of particle current, with w(,(ft 7 ) playing the role of 
the (relative) group velocity of the asymptotic mode 7 . 
As any linear combination of eigenstates with the same 
energy oj is again an eigenstate, this relation can be ex¬ 
tended to 


This “effective scattering matrix” consists of the first 
(M + N) x (M + N) blocks of the exact effective Hamilto¬ 
nian and norm matrix and the A matrices of the asymp¬ 
totic part [Eq. ©] to make sure that these matrices 
remain the truncated versions of a hermitian problem. 
This matrix has D 2 {d — 1) + (N + M)L rows, which 
implies that the linear equation ( |17[ ) has T exact solu¬ 
tions, which is precisely the number of scattering states 
we expect to find. Every solution consists of a local part 
( D 2 (d — 1) + NL elements), the LM — T coefficients p 
of the decaying modes and the 2 T coefficients q of the 
asymptotic modes. 


F. S matrix and normalization 

After having shown how to find the solutions of the 
scattering problem, we can now elaborate on the struc¬ 
ture of the asymptotic wave function and define the S 
matrix. 


Y,<l]( K ,u)qU K ^) (~^( K 'i')) =0 - 

With this equation satisfied, we can define the two- 
particle S matrix S(K,co). Firstly, the different modes 
are classified according to the sign of the derivative: the 
incoming modes have 77 > 0 (two particles moving to¬ 
wards each other), the outgoing modes have 77 < 0 (two 
particles moving away from each other), so that we have 


E q]{K,u)q7(K,u,) 

7er in 


dw 7 

dre 



= E <i]( k ,^)qUK,uj) 

7er out 


dw 7 

An 



If we group the coefficients of all solutions in (square) 
matrices Q m {K, u) and Q 0 ut{K, w), so that the i’th col¬ 
umn is a vector with the coefficients q] for the in- and 
















outgoing modes of the i’th solution, we can rewrite this 
equation as 


Qin(K, u)^(K, w)Qin(K, u) 

= Qout(K, w) t Ct {K, to)Q out (K, u ). 


with Vin.out (K,u>)ij = 5ij ^A(k 7 ) a diagonal matrix. 

As Qi n (K,cj) and Q on t{K,uj) should be connected lin¬ 
early, we can define a unitary matrix S(K,lj) as 


1/2 


H out (A»Q out (A» = SiK^V^K^Q-^K^). 

In Sec.|III A|we will show that this definition corresponds 


to the standard S matrix. Note, however, that S(K,ui) 
is only defined up to a set of phases. Indeed, since the 
vectors v 1 can only be determined up to a phase, the 
coefficient matrices C- In and C 0 ut are only defined up to a 
diagonal matrix of phase factors. These arbitrary phase 
factors show up in the S matrix as well. We will show 
how to fix them in the case of the elastic scattering of 
two identical particles (Sec. |HG| ; in the case where we 
have different outgoing channels only the square of the 
magnitude of the S matrix elements is physically well- 
defined (see Sec. Ill A l. 


This formalism allows to calculate the norm of the scat¬ 
tering states in an easy way. Indeed, the general overlap 
between two scattering states is given by 


(TAK^oj'^T^K^)) =2tt5(K - K') ( ^ q^'(K f , w')?7(tf, u)v\,v y ^ e i( ^“V )n + finite 

17 , 7 ' n,n'>iV 


= 27 t6(K - K') q?' (K',u')q?(K, u)v\,v^7rS(K^(u;) - «y(u/)) + finite . 




The S factor for the momenta k 7 is obviously only satisfied if u> = a/, so we can transform this to a 8(ui — u>'). 
Moreover, if k 7 (w) = k! for 7 7 ^ 7 ', then necessarily VyV 7 = 0, so we can reduce the double sum in 7 , 7 ' to a 
single one. If we omit all finite parts, we have 

(T»/ (K 1 ,uj')\Ti(K,uj)) = 2-ir5(K - K')ttS(uj - w') ^ ^(K, w) 


dw 7 

d« 


(k 7 


With the Qi n /out as defined above the overlap reduces to 

(Y.j' (K 1 , oj)) = 2n5(K - K')2nS(uj - J) [Q in (K,V?(K,w) [Q in (K,uj)] l 

= 2nS(K - K')2tt8(uj - w') [Q out {K,Lj)}\, V 2 ut (K,w) [Q„ut (!*»]< - 


G. One type of particle 

Let us make things more concrete by working out the 
case where the one-particle spectrum consists of just one 
type of particle with dispersion relation A(k). Suppose 
we have only one combination of momenta K\ and k 2 
such that they add up to total momentum K = «q + k 2 
and total energy oj = A(ki) + A (k 2 ). There are two 
asymptotic modes - one mode with kj on the left and 
K 2 on the right, and one mode with the momenta inter¬ 
changed - that combine into one scattering state with 
the asymptotic form 

c(n) q 1 e iK2n v 1 + q 2 e iKin v , 2 . 

The conservation equation that was derived in the previ¬ 
ous section takes on the simple form 



because u/(k 1 ) = —u>'(k 2 ). As we mentioned above in 
the general case, the relative phase of the two vectors iq 


and V 2 can be chosen arbitrarily. However, since the two 
modes correspond to the interchanging of two identical 
particles, it makes sense to fix the phase such that v\v^ > 
0. Due to the momentum dependence of the one-particle 
solutions, this overlap will be slightly smaller than one. 

The S matrix reduces to a phase factor and is defined 
as 

S(k i,k 2 ) = S(K,ui ) = 

T 

The asymptotic wave function takes the form 

|T(ki,k 2 )) -> ^2 e * (Kini+K2 " 2) [B Kl at ni,B K2 at n 2 ] 
" 1<2 

+ S(« 1 , K 2 )e i(,t2ni+Kin2) [ B K2 at ni ,B Kl at n 2 ]. (20) 

From simple arguments'^ one can argue that in one di¬ 
mension the S matrix should have t he universal limiting 
value for low-energy scattering® 7 ! 

S(ki, k 2 ) —► —1 as |/«i — k 2 \ —> 0. 
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We define the scattering phase 9 as the phase shift of 
the S matrix relative to its universal low-energy value 
S(m,K 2 ) = ~e ie ^’ K2 \ 

H. Spectral functions 

With the variational wave functions of one- and two- 
particle states, we can now calculate the low-energy part 
of spectral functions at zero temperature. We consider 
the following function 


and 

«i(k) = J ^ <T 0 | Ol K 2n6(u) - H)O 0 ( 0) |*„> 

= <*o|o! KJ ffOo( 0 )|tf 0 >. 

If the ground state is taken to be an MPS, these quan¬ 
tities can be calculated exactly. Note that so is just the 
static correlation function and that the ratio of the two is 
equal to the single mode approximation for the dispersion 
relation^- 


S{k,u>) = [ d<tf 0 | Ot (i)O 0 (0) |* 0 > 

n 

with O n (t) an operator at site n in the Heisenberg pic¬ 
ture. In order to approximate the low-energy part, we 
insert a projector on the one- and two-particle subspaces 



where Id (r 2 ) is the set of all types of one-particle (two- 
particle) states at that momentum (momentum-energy). 
The states are orthonormalized as 

(* 7 /(«/)|* 7 (k)) = 27 tS(k' — k)5 77 ' 

(T 7 / (K 1 , u/)|T 7 (A', ui)) = 4tt 2 (5(A" - K)S(u’-u)S rf 

so that we obtain the Lehmann representation®! for the 
spectral function up to two-particle contributions 


S(k,u)= ^2 27t5(A C[ (k) - w) ($£*(«)| Oo |*o) 

aGT^/t) 

+ ^2 w )l 1 * 0 / 

+ ... 


In gapped systems, the one- and two-particle contri¬ 
butions saturate the full spectral function below the 
three-particle threshold, while contributions from higher- 
particle excitations might become important for higher 
energies. Yet, it appears that typically the one- and two- 
particle sectors already contain the largest portion of the 
spectral function, see e.g. Ref. H§J The one- and two- 
particle form factors appearing in the above expression 
are calculated explicitly in Appendix |A 3| 

To get a quantitative estimate of how well the spectral 
function is approximated, we look at the zeroth and first 
frequency moment at a certain momentum, which are 
defined as 


so(k) = J —S(k,uj) and Si(k) = J —u>S(k,u>). 

These quantities follow the sum rule^l 


«o(«) = J ^ <*o| Ol K 2n6(oj - H)O 0 (0) |T 0 ) 

= <*o|Ol*Oo(0)|*o> 


\ , , Sl ( K ) (*o\Ol K H0 0 (0)\*o) 

Asma(k) - ao («) " <* 0 |Ol B Ob(0)|*o) ' 

By comparing the one- and two-particle contributions for 
so and Si to the exact values, we can get an idea of how 
well these eigenstates capture the effect of the operators 
working on the ground state and, consequently, how well 
the spectral function is approximated by only looking at 
these contributions. 


III. TWO-PARTICLE S MATRIX AND 
APPROXIMATE BETHE ANSATZ 

We now discuss how the variational formulation of 
scattering theory using matrix product states, as devel¬ 
oped in the previous section, relates to standard scat¬ 
tering theory. We then discuss how we can use the in¬ 
formation provided by the scattering matrix to build an 
effective description of the low-energy behaviour of the 
spin chain using the approximate Bethe ansatz. 


A. Stationary scattering states and the S matrix in 
one dimension 

In standard scattering theory the S matrix is typically 
defined from a dynamical point of view: its elements are 
the overlaps of asymptotically free in and out states with 
respect to the full time-evolution operator. Although it 
is a priori not clear that this definition corresponds to 
the one that was presented in the previous sections, we 
can show that this is indeed the case. 

Appendix [D] provides a summary of the standard scat¬ 
tering formalism of single-particle quantum mechanician, 
which we have adapted to the one-dimensional setting 
with general Hamiltonians (e.g. potentials which are not 
diagonal in real space) and arbitrary dispersion relations 
(non-quadratic eigenvalue spectrum of the “free” Hamil¬ 
tonian Hq). More specifically we have shown how the 
S matrix elements f(qp <— p a ) show up in the asymp¬ 
totic form of the scattering eigenstates |p a ±) of the full 
Hamiltonian H. 

To make the connection to the variational scattering 
states of Sec. [TTJ we have to make a few modifications. 
First of all, we can reformulate the two-particle scattering 
problem as a one-particle problem by factoring out the 
conservation of total momentum and only focus on the 
matrix elements between different relative momenta. At 
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every value of the total momentum, we can define relative 
momentum states |p 7 ) with dispersions u/(p 7 ), which are 


solutions of the free Hamiltonian Hq. This free Hamil¬ 
tonian corresponds to the effective two-particle Hamilto¬ 
nian matrix in the asymptotic regime (111 and the states 
| pj) are the asymptotic modes (14). 

Secondly, our “one-particle” Hilbert space is only de¬ 
fined on a half-infinite line, because the particles are es¬ 
sentially bosonic. The way around this consists of ar¬ 
tificially assigning particle labels and distinguishing the 
situation where particle 1 ( 2 ) is on the left (right), and 
the opposite situation; the relative coordinate n = n 2 —n\ 
now ranges over the positive and negative integers. Al¬ 
ternatively, one could add to the free Hamiltonian Hq 
a potential V which is infinite everywhere on the nega¬ 
tive real line, making this a forbidden region. Scattering 
theory would still work (existence of the Moller opera¬ 
tors etc.), provided that we restrict the “in” states to 
momenta for which ^(p) < 0 and the out states to mo¬ 
menta for which ^(p) > 0. This corresponds exactly 
to how we defined the incoming and outgoing modes in 
Sec. ITTFl 

Translating the expression for the asymptotic wave 
function of the scattering states |p a +) to the framework 
of Sec. m amounts to the following form for the wave 
function c a (n ) 


c a (n) 


£(*.) 

dft 


- 1/2 


v a e lp “ n 


Y /(« 7 K a ) 

jGA+(K a ) 




- 1/2 


v~e 


for every incoming mode a = 1,..., T. In this represen¬ 
tation, we choose one incoming mode a that couples only 
to all outgoing modes {7 £ A + (/c a )}. The coefficient ma¬ 
trix for the incoming modes that was defined earlier takes 
on the form 

- 1/2 

{Qin)j,ct — ^70 




while the coefficients for the outgoing modes are given by 


(Qout)y,a 


dw 

d/c 


(/c 7 ) 


- 1/2 


f (/c 7 ^ /c a ). 


The S matrix S(K,co) that was defined takes on the sim¬ 
ple form (as Vi n Q m = 1 in this representation) 

S{K,u) = V out Q out 

= /(« 7 < r - K a ). 

Through this identification the unitariness of the S ma¬ 
trix S ( K , ui) that was proven in the previous section is 
indeed equivalent to the unitary S matrix defined through 
the Mpller operators as S = 


B. Scattering length and bound states 

Suppose we have the scattering process of two identi¬ 
cal particles in the limit of vanishing relative momentum. 


We expect that the equation for the relative wave func¬ 
tion ij){x) should obey the zero energy and zero potential 
Schrodinger equation 


d 2 4’{x) = 
dx 2 

in the region x > Xq where Xq is the length of the inter¬ 
action. The solutions are of the form oc x — a for 
large x which matches the asymptotic form of Sec. |II G| 
if the phase of the S matrix reduces to 

0(«4, «2) ~ —a(«i — K 2 ) (21) 

in the limit for ki — n 2 —> 0. The slope a will be called 
the scattering length and still depends on the total mo¬ 
mentum Hi + K 2 - 

Suppose now the existence of a bound state with very 
low binding energy — e. The wave function of this bound 
state should look like i/j(x) = e~ KX ~ 1 — kx with 
ui(in) = — e —> 0. If we want the formation of this 
bound state to follow smoothly from a scattering state 
with vanishing energy, the scattering length should di¬ 
verge. This means that the formation of a bound state 
out of a scattering continuum at a certain momentum 
should be accompanied by a diverging scattering length. 


C. Approximate Bethe Ansatz 

In this section we will develop a method to describe 
a finite density of excitations based on the coordinate 
Bethe ansatz. For simplicity, we will for the remainder 
of this section restrict to the case of one type of particle - 
making the consistency conditions for factorized scatter¬ 
ing (Yang-Baxter equation) trivial - but the framework 
can be extended to multicomponent situatiomPSI. We will 
interpret the strongly correlated MPS ground state as a 
vacuum state on which we can build iV-particle states, de¬ 
scribed by a TV-particle wave function \H(a:i,..., Xn)- Al¬ 
though in general we have no particle conservation in the 
system, we will argue that the first-quantized approach 
gives a good approximation at low densities. Indeed, 
particle-number violating processes involve three or more 
particles and can be neglected at low densities. In Sec. [V] 
we will discuss how to develop a second-quantization ap¬ 
proach. 

We start with one particle. We can link the one- 
particle excitation |$ K [5]) with dispersion A (k) in an 
obvious way with a one-particle wave function 4'i(a:) in 
first quantization as 


4b (aO =e iKX . 

Adding a second particle can be done by only taking 
account of the asymptotic part of the two-particle wave 
function [Eq. p0| )] (x\ < X 2 ) 

^ 2 (x 1 ,x 2 ) = e^ KlXl+K2X ^ + S{ Kl , K 2 )e^ K2Xl+K1X2 \ 

As we are working with identical particles, the wave func¬ 
tion in the other sector (aq > x 2 ) has to be determined 
from the statistics of the particles. On the level of the 
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spin system, the addition of a particle is a local opera¬ 
tion, so we will work with bosonic many-particle wave 
functions. 

The addition of a third particle can only be done ap¬ 
proximately. Indeed, a three-particle wave function has 
the general forirPS 

* 3 (*i, *2,0:3) = e i(KlXl+li2X2+K3X3) 

+ K 2 )e i{K2Xl+KlX2+K3X3) 


d«id«4d*& 5(444) e i ^ Xl+ ^ X2+ < X3) 


other particle numbers. 


( 22 ) 


The first terms represent a sum over all six permutations 
of the three momenta, with the S matrices for all possible 
two-particle scattering processes as prefactors. The next 
term is the diffractive part, which accounts for the three- 
particle scattering. For these scattering processes, the 
two conservation laws are not enough to preserve indi¬ 
vidual momenta and we can generate a whole continuum 
of other momenta. The last term accounts for the non¬ 
particle preserving scattering processes, which can gen¬ 
erate two- or four-particle states as well. As a result, it is 
no longer possible to assign a set of individual momenta 
{«i, « 2 , K 3 } (or even a particle number) to this wave func¬ 
tion, because they are completely mixed up with all other 
possible sets of momenta that are compatible with con¬ 
servation of total energy and momentum. 

The crucial approximation of our approach is that we 
neglect the two last terms in Eq. (22): every many- 


particle scattering event can be decomposed into two- 
particle scatterings that preserve particle number and 
individual momenta. This implies that three-particle 
eigenstates can be labeled by three individual momenta 
and that the three-particle wave function is given by 
the permutation terms only. The absence of diffrac¬ 
tive scattering is the hallmark of integrabilitjEl, so we 
are essentially assuming that our many-particle system 

is integrabld^HIZl 

If this approximation proves to be valid, we can ap¬ 
ply the Bethe ansatz machinerjESMU. The first-quantized 
wave function of an integrable IV-particle system, unam¬ 
biguously defined by a set of momenta {Ai,..., Am}, is 
a sum of plane waves with all possible permutations of 
the momenta 

^(x u ...,x N ) = J2A{P)e i ^' plXl+ ' +x ' PNXN) (23) 
v 

where A(V)/A(V') = 5(A,, A j) if the permutations V and 
V' differ by the interchange of the momenta Aj and A j. 

By imposing periodic boundary conditions on the 
Bethe wave function in the thermodynamic limit, we ar¬ 
rive at a description of the ground state as a Fermi sea 
of “pseudo-momenta” filled up to a certain Fermi level 
q. In contrast to the free-fermion case, the density of oc¬ 
cupied modes is not constant but given by the function 
p{ A) such that p{ A) = 0 for |A| > q. The energy of the 
modes e(A) can be determined from the integral equation 


e(A)- 


1 

27T 


where eo(A) is the “bare energy” of the particle, i.e. the 
energy an isolated particle with momentum A would have 
in an infinite system. The kernel of the integral equa¬ 
tion is given by the derivative of the scattering phase 
K(X 1 p) = d\0{X,p). The value of the Fermi level is 
computed self-consistently from this equation and the 
requirement that e(±g) = 0. Once q has been deter¬ 
mined, the density p( A) is the solution of a similar inte¬ 
gral equation— 

p ( A ' ) — 2 bj = 24 ( 25 ) 

The total density and energy density are given by 

D= r p(A)dA and E = — T e(A)dA. (26) 
J-q 27T J_ q 

The excitation spectrum is easily characterized in 
terms of the pseudo-particles of the Bethe ansatz. We 
can construct two types of elementary excitations: ei¬ 
ther we take one particle with momentum |A| < q out 
of the Fermi sea (hole excitation) or we add one particle 
with momentum |A| > q (particle excitation). These ele¬ 
mentary particle and hole excitations have a topological 
nature 22, so that the physical excitations - the ones hav¬ 
ing a finite overlap with a local operator - consist of an 
even number of particles and liolesP^. This gives rise to 
the physical excitation spectrum as shown in Fig. |5(b)| 

This critical one-dimensional bose gas can be described 
as a Luttinger liquid (LL)P25 7 . A first important quan¬ 
tity is the Fermi momentum kp, the physical momentum 
of the gapless particle and hole excitations. It is given by 
the dressed momentum of the Fermi level and is directly 
related to the density as (see appendix) 


Uf = 7T D. 


(27) 


Since we have gapless excitations at 0 and A2kp. cor¬ 
relation functions will have their oscillation periods at 
these values. The slope of the dispersion relation is the 
Fermi velocity u and can be calculated from the Bethe 
ansatz. The third important characteristic quantity is 
the LL parameter I\ which determines the power-law de¬ 
cay of correlation functions. In order to calculate it, we 
define the function S r(A) as (h is a chemical potential 
for the particles) 


Sr( A)= - 


fMA) 

dh 


which (from Eq. (24)) follows the integral equation 


5 fl (A)-^ j" K(X,p)S R (p)dp 


= 1 . 


(28) 


K( A, p)e(p)dp = e 0 (A) (24) 


In the context of a dilute gas of magnons (see Sec. IV D) 
S R (q) can be interpreted as the renormalized spin of the 
magnon close to the Fermi surface. With the low-energy 
excitations just above the Fermi sea behaving as free 
fernhoniP^ (i.e. their S matrix is -1), one can show that 
the LL parameter I\ is related to S R (q) aJ^l 

K = S R (q) 2 . 


-q 


(29) 
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integral equation zero, we find easily the density and the 
(physical) Fermi momentum 



and the LL parameters 

u = 27 tcD, K = 1. 


Upon increasing the density, the limiting value of the 
S matrix will no longer apply. From Sec. eh B| we know 
that the first order correction to the scattering phase 
is given by the scattering length, so we can insert the 
form © into the Bethe equations, while still assuming 
a quadratic dispersion relation. The first order correction 
to the Fermi level is linear in the scattering length 


Q = Qff + Sq = q ff — 


ah 
3i re’ 


(31) 


so that the correction to the density is given by 


D = 


1 


4 ha 
37r 2 c 


+ 0(a 2 ). 


(32) 


FIG. 5. (a) The Fermi sea of pseudo momenta, filled up to the 
Fermi level q. Physical excitations can be pictured as particle- 
hole excitations close to the Fermi-level, (b) The physical ex¬ 
citation spectrum, the grey area represents a continuum of 
states. Because of the fact that physical excitations always 
come in pairs, the spectrum has its minima at momentum 0 
and 2k,F- The slope of the dispersion relation at these mo¬ 
menta is the Fermi velocity u. 


This result coincides with the one in Re f. (521 The LL 
parameters in first order in a are given bj^ 63 ^31 

„ h 4ah _. 9 , 

u = 2 c\j —I——-b 0(a ) 

V c 37 t 

and 

K = 1-2 aD + 0(a 2 ). 


By thus making the connection between the approximate 
Bethe ansatz and the LL description, we can infer infor¬ 
mation on the critical correlations in a system where a 
finite density of excitations forms on top of a strongly- 
correlated vacuum state. More specifically, we can infer 
the long-range behaviour of one-particle and pair corre¬ 
lation functions ad—! —! 


9i(x) = A 0 


1 

X 1 /2K 


— Ai 


cos(2ttDx) 

X 2K+1/2K ~*~ 


D 2 {x) = D 2 


K cos(2nDx) 

2n^ + 1 X™ + 


(30) 


where D is the density, A 0 , A\, and Bi are non-universal 
constants and the dots denote higher order terms. De¬ 
pending on whether the operator targets a particle or a 
pair, the corresponding correlation functions will decay 
according to one of these two forms. 


D. Limiting cases 

The Bethe ansatz equations of the previous section 
can be greatly simplified if we assume that we work at 
very low densities. Indeed, assuming that only the lowest 
pseudo-momentum states are occupied, we can approx¬ 
imate the full dispersion relation by its quadratic form 
eo(A) w cA 2 — h, and the full two-particle S matrix by its 
limiting value of S(6,p) ss —1. With the kernel of the 


E. Thermodynamic Bethe ansatz 


At zero temperature, the coordinate Bethe ansatz de¬ 
scribes an integrable system in its ground state by filling 
up a Fermi sea of quasi-momentum states; its excitations 
are holes and particles above this Fermi sea. When a 
finite temperature T is applied, these particles and holes 
will have finite distribution densities. By associating an 
entropy to these distributions and minimizing the free en¬ 
ergy, one arrives at the celebrated Yang-Yang equation — 1 


e(A) = e 0 (A) 


T 

2n 


[ + K( A,^)ln(l + e- e( ' i)/T )d/r, 

J —OO ^ ' 


a non-linear integral equation for the dressed energy e(A) 
of the quasi momentum states; the equation can be solved 
by iteration®!. The density of occupied vacancies p( A) is 
given by 


0{ A) 


P{ A) 

Pv( A) 


1 

\+e< x )/ T 


with p v ( A) the density of all (occupied and empty) va¬ 
cancies. Through this equation the density of occupied 
vacancies satisfies the integral equation 

p(A) = ^^(i + / K(X,p,)p(p)dpj , (33) 
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such that the total density can be calculated as 


D = 



p{ A)dA. 


(34) 


F. Effective integrable field theories 



Another way of dealing with a finite density of excita¬ 
tions, based on information on the one-particle dispersion 
and the two-particle S matrix, consists of mapping the 
system to an effective integrable field theory. The pa¬ 
rameters in this effective theory should be tuned to fit 
the variational information as good as possible. This ap¬ 
proach has the advantage that integrability is exact for 
the effective model, but the mapping is typically only 
valid in some small region (e.g. low density and/or low 
temperature). 

One possible field theory is obtained by making the 
approximation that the particles interact through a con¬ 
tact potentia l 37 ! 67 !, so that we end up with a Lieb-Liniger 
modeP^. The first-quantized Hamiltonian for a collection 
of N bosons is given by 



with the mass m of the bosons and the interaction 
strength c as the two tunable parameters. The two-boson 
S matrix is given by <S(Ai, A 2 ) = — e *0(Ai-A 2 ) w ith 


9( A) = 2 arctan 



(36) 


so that the scattering length for a 6 potential is as = 
— 2/c. The boson dispersion relation is just quadratic, 
i.e. A(A) = A 2 /(2 m). By variationally calculating the 
dispersion relation and the scattering length of the rele¬ 
vant excitations, we can fix the two parameters and map 
the density of excitations to a Lieb-Liniger model. At low 
densities, we expect that this mapping is quantitatively 
correct. 

Another possibility is the non-linear sigma model, 
which has proven to capture the qualitative behaviour 
of Haldane-gapped spin chains such as the spin-1 Heisen¬ 
berg modeP^or two-leg spin- 1/2 ladders. In contrast to 
the Lieb-Liniger model, however, we can not tune any 
parameters to fit the exactly known 6 ® two-particle S ma¬ 
trix. The universal behaviour of e.g. the magnon con¬ 
densation of a gapped spin chain in a magnetic field®!, 
can nonetheless be captured with this mapping. 


FIG. 6 . The ladder geometry with Jy and Jj_ the couplings 
along the leg, resp. rung. We will always put Jm = 1 and 
define the coupling ratio 7 = 


where l = 1,2 denote the two legs of the ladder and 5/.; 
denotes the spin operator at site i in the Z’th leg (see 

Fig.| 6 |. 

The two-leg HAF ladder and its excitation spectrum 
have been studied intensively for many reasons. First 
of all, it is the first step to study the transition from 
one-dimensional systems to higher-dimensional versions. 
Secondly, the excitation spectrum has a lot of interest¬ 
ing features, such as the presence of a gajP 6 and the 
existence of bound states, and can be studied with a 
variety of methods depending on the parameter regime. 
These features can also be observed experimentally' hEU, 
so that lad ders p rovide an ideal test for these theoreti¬ 
cal methods ' 5 76 4 Finally, the experimental realization of 
magnetized spin ladders provides an ideal quantitative 
test of the Luttinger liquid model 77 79 . 

In this section we will test our variational method 
on the two-leg ladder. An MPS approximation for the 
ground state can be found by first blocking two spins 
on a rung into one four-level system and applying an 
MPS opt imiza tion algorithm (we have used the TDVP 
algorithm ! 32 ! 33 !). In this representation (to every rung 
there corresponds one MPS tensor A) we find a ground 
state that is invariant under translations over one site in 
the leg direction; all momenta in the following subsec¬ 
tions are defined with respect to this translation opera¬ 
tor. The Hamiltonian and the ground state are invariant 
under the reflection operator V which flips the two legs 
of the ladder. We impose no additional symmetries (e.g. 
SU(2) invariance) on the MPS, but our variational so¬ 
lution will of course have the right symmetries to high 
precision. 

In the first three subsections we will investigate the 
low-lying spectrum of the ladder without magnetic field. 
In the following two subsections we will apply the ap¬ 
proximate Bethe ansatz to the magnetization process, at 
zero and finite temperature. 


IV. APPLICATION TO SPIN LADDERS 


A. One-particle excitations: elementary spectrum 
and bound states 


We will study the spin-1/2 Heisenberg antiferromag¬ 
netic (HAF) two-leg ladder in a magnetic field, defined 
by the Hamiltonian 

H = J2 1 ■ Si + i,i +7 X] ^-1 • Si,2 Sh (37) 

i,l i i,l 


The nature of the elementary excitations in the ladder 
can be understood starting from different limits. 

At zero coupling (7 -A 0), we have two independent 
spin-1/2 Heisenberg chains where the elementary excita¬ 
tions are spinons (carrying spin 1/2). These spinons are 
topologically non-trivial excitations and can only be ere- 
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momentum 


FIG. 7. The rescaled gap ^ A 2 in function of the interchain 

coupling 7 . The blue dashed lines are the first order correction 
from the strong-coupling limit (7 —> 00 ) and resu lts from 
bosonization in the weak-coupling limit (7 —> q^ 80 I 83 [ 


ated in pairs by the action of a local operator. Upon cou¬ 
pling the chains, the spinons are confined into magnons 
carrying integer spin. This picture has been studied with 
bosonization technique^} showing that the interchain 
coupling opens up a gap to a triplet of massive magnons 
(triplons) and a higher-up singlet. 

At infinite coupling (7 —► 00 ) we have a collection of 
independent rungs with antiferromagnetic interaction. In 
the ground state all rungs are in a singlet state and an 
elementary excitation is constructed by promoting one 
rung to a triplet state. When the leg coupling J± is 
turned on, this triplet obtains a kinetic energy and we 
get a non-trivial dispersion. This qualitative picture sur¬ 
vives for intermediate couplings: through perturbative 
continuous unitary transformations an effective particle 
picture can be established and very accurate results on 
e.g. the elemen tary d ispersion relation and bound states 
can be obtained ^ — ' 1 

In Fig. [7] we have plotted the gap in function of the 
interchain coupling. One can observe that the gap goes 
to zero in the weak-coupling limit, while it grows to the 
constant value that one expects from a strong-coupling 
expansion. Our variational results smoothly interpolates 
between these two limits. 

A typical excitation spectrum in the intermediate re¬ 
gion (7 = 2) is shown in Fig. [ 8 ] The lowest-energy state 
is an elementary triplet excitation (magnon) with a min¬ 
imum at momentum n. The magnon has odd parity un¬ 
der the reflection operator V. The lowest-energy state 
around momentum zero is a two-magnon scattering state 
and has even parity. Because the one- and two-magnon 
state have different parity, the elementary magnon can¬ 
not decay and is stable in the whole Brillouin zone. From 
Fig-! where we have plotted the variance of the excita¬ 
tion ansatz, we can indeed see that the magnon is a bona 
fide particle excitation for all momenta. Note that under 
a parity-breaking interaction the stability of the magnon 
inside the continuum breaks down 84 and it might prove 
an interesting question whether we can capture its decay 
within our framework. 


FIG. 8 . The one-particle spectrum consists of a triplet 
(magnon) which is stable over the whole Brillouin zone (lowest 
lying blue curve), an singlet (bound state) which is stable for 
momenta between kbsi ~ 0.397T and 7r (second blue curve), 
and a triplet (bound state) which is stable for momenta be¬ 
tween kbs 2 ~ 0.46-7T and n (third blue curve). Note that the 
determination of kbsi and kbs 2 is not very precise because 
the one-particle ansatz is not accurate near the transition. 
The red region is the two-magnon continuum, the green re¬ 
gion is the three-magnon continuum; the other continua (e.g. 
triplet-singlet continuum) are not shown. 


The elementary excitation spectrum at 7 = 2 has 
two more elementary particle excitations, a singlet and a 
triplet, which are stable in a limited region around mo¬ 
mentum 7 r. Both are even under the parity operator 
V. From the strong-coupling expansion, we can inter¬ 
pret them as two-magnon bound stated®} hence the even 
parity (without a well-defined particle number, we can¬ 
not make this interpretation, so we regard these branches 
as elementary particles). The variance of the bound 
states is sufficiently small in the stable region, but it 
grows larger as the momentum approaches the contin¬ 
uum. From Ref. [5] we know that the localized nature 
of an elementary excitation is related to the gap below 
and above the excitation branch, so we expect the bound 
state to become wider as the gap to the continuum closes. 
This explains the increasing variance of the bound states 
in Fig. [9} Upon entering the continuum, the bound state 
has become completely delocalized and no longer exists 
as a stationary eigenstate of the Hamiltonian. 

As a last illustration of the one-particle ansatz we have 
included Table[l]with excitation energies and variances in 
the weak-coupling region, showing the elementary triplet 
and singlet excitations that we expect from a bosoniza¬ 
tion calculation. We observe that the variances are some 
orders of magnitude larger in this weak-coupling region. 
Since the gaps above and below these excitations are a 
lot smaller at small 7 , this is not unexpected. Note that 
both the energies and the variances have the right degen¬ 
eracies, even though we never imposed the corresponding 
symmetries explicitly. 
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FIG. 9. The (log 10 of the modulus of the) variance of the 
one-particle excitations; dots, resp. crosses are positive, resp. 
negative variances (see Appendix |A 4| for the meaning of a 
negative variance). The magnon (green) is clearly a well- 
defined particle excitation in the whole Brillouin zone. The 
singlet (red) and triplet (blue) get larger variances as they 
come closer to the two-particle band, until they actually dive 
in and are no longer stable. Calculations were done at 7 = 
2 with bond dimension D = 30; the ground state variance 
density is 2.27 x 10 -8 at that bond dimension. 


energy 

variance 

0.081841224772803 

-0.000178252361115 

0.081841224779434 

-0.000178252351941 

0.081841224792513 

-0.000178252347304 

0.331378942771407 

0.000337897356458 

0.367322866763615 

0.029803975299627 

0.410460620351393 

0.044970779553592 



0.513408977989184 

0.014052233372105 

0.513408978649963 

0.014052233100514 

0.513408978939573 

0.014052232922150 




TABLE I. Excitation energy and variance of the first 6 solu¬ 
tions of the one-particle problem for the HAF (7 = 0.2) at mo¬ 
mentum 7 r with bond dimension D = 108. The variance den¬ 
sity of the ground state is 9.28.10 -6 . The first triplet has neg¬ 
ative variance, which shows that this excitation is closer to an 
exact eigenstate locally than the ground state (see Appendix 


A 4 1 . The fourth solution is also a true one-particle (singlet) 


excitation. All other solutions have a considerably larger vari¬ 
ance and correspond to artificial two-particle states. Further 
up in the continuum, however, we have another triplet with 
quite small variance, although it is difficult to say whether 
this corresponds to a true bound state. 


B. Two-particle S matrix 


In this section we will look at the two-magnon S ma¬ 
trix; the scattering of, e.g., an elementary magnon with a 
bound state will not be considered. The S matrix was de¬ 
fined in Secs. |TTTland |mAl in our setting we have three 
types of particles (the three components of the magnon 
triplet) and they all have the same dispersion relation. 
This implies that, for every combination of total momen- 


FIG. 10. The S matrix in function of relative momentum 
Ki — K 2 at total momentum K = 0. Plotted are the phases 
of the S matrix in the S = 0 (red), S = 1 (blue) and S — 2 
(green) sector. Calculations were done at 7 = 2 and with 
bond dimension D = 32. 


turn K and total energy w within the two-magnon con¬ 
tinuum, we can build 9 scattering states. The relative 
coefficients of the asymptotic modes in these scattering 
states give rise to a (9 x 9) unitary S matrix (the group 
velocities will factor out, as all particles have the same 
dispersion). Furthermore, instead of labeling these scat¬ 
tering states with momentum and energy (K, w), we can 
equally well label them with total and relative momen¬ 
tum (K, K\ — K 2 ) where K\ and K 2 are the two momenta 
that show up in the asymptotic modes (there is still an 
ambiguity in the ordering of the momenta, we will always 
take the convention that k± > « 2 , he. positive relative 
momentum). 

We can simplify the S matrix by making use of SU(2) 
invariance. Indeed, if we make linear combinations of 
the asymptotic modes that diagonalize the total spin 
and its projection S f., the S matrix should be diagonal. 
Moreover, since the magnon interactions are SU(2) in¬ 
variant (both Hamiltonian and ground state are), the S 
matrix elements should be constant within every sector 
of total spin. This means that the general expression 
for the magnon-magnon S matrix in this representation 
should reduce to 

/ e * e ° x l lxl 

s = ~e iSl x l 3x3 

V —e 102 X 15x5 



i.e. the S matrix reduces to three phases for every sector 
of total spin. In our simulations, we always found this 
reduced form to high precision, so in the following we can 
restrict to plotting these three phases. 

In Fig. [TO] we have plotted the S matrix in function 
of the relative momentum ki — K 2 for total momentum 
K = 0. One can observe (i) the limit S = — 1 for the rel¬ 
ative momentum going to zero, and (ii) the linear region 
around this limit (the slope is the scattering length). The 
sign of the phase is positive for all three sectors (although 
this does not have to be the case, see Figs. 11 and[l2|). 

In Fig. El we have plotted the S matrix in the S = 2 
sector for different values of the total momentum. We 
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relative momentum 

FIG. 11. The scattering phase in the S = 2 sector for 8 
equally spaced values of the total momentum between K = 
0 (upper line) and K = 7t/3 (lower line). Around K = 0 
there is a region where the S matrix is independent of total 
momentum, which points to some Galilean invariance around 
the minimum of the dispersion relation. Calculations were 
done at 7 = 2 and with bond dimension D = 32 



FIG. 12. The scattering lengths ao (red), an (blue) and 02 
(green) in function of the total momentum K. In the S = 2 
sector nothing spectacular happens, although it does change 
sign. In the other sectors we see a divergence at the momen¬ 
tum where a bound state forms. The plotted range does not 
show all data points around the divergences, the full lines are 
a guide to the eye and give an indication on where the other 
points are situated. Calculations were done at 7 = 2 and 
bond dimension D = 32. 


observe that the S matrix depends strongly on K in a 
non-trivial way, but there seems to be a small region 
around K = 0 where it is quasi-constant. This points 
to the presence of a region around the minimum of the 
dispersion relation where the interaction is Galilean in¬ 
variant (note that the dispersion should be quadratic in 
this region). At larger momenta, this Galilean invariance 
is broken, as one expects in a lattice system. 

Even more spectacular things can happen when we 
vary the total momentum, such as the formation of a 
bound state. In Fig. [12] we have plotted the scattering 
lengths in all three sectors in function of the total mo¬ 
mentum. We can see that the scattering lengths in the 
S = 0 and S = 1 sectors diverge, signalling the formation 
of the singlet and triplet bound states (in agreement with 
the discussion in Sec. Ill B). 


C. Spectral function 

Since we have a two-leg ladder system, we can look at 
spectral functions with transversal momentum q equal 
to 0 or 7 r. We define the two rung operators (defined on 
rung *) 


(SS)i = Sl 1 + St i2 (38) 

(S~) t = SI 1 - S? 2 . (39) 

These operators have even, resp. odd parity under the ac¬ 
tion of the reflection operator V. We will look at spectral 
functions So/ n (n, w) with respect to these two operators, 

<V(k,w) = £ / die !( “‘- Kn) 

n J 

x <*o| e- iHt {S* 0 /J n e iH \S * 0/n ) 0 |tt 0 > (40) 

where ^2 represents a sum over rungs. 


Let us first look at the one-particle contributions. 
Since the elementary magnon is odd under V, it can only 
carry spectral weight with respect to the odd operator. 
From SU(2) symmetry we know that the singlet bound 
state does not carry any spectral weight with respect to 
both operators (they are both spin-1 operators). Lastly, 
the triplet bound state is even under V, so it only con¬ 
tributes to the even operator spectral function Sq(k,ui). 
These considerations lead to the picture in Fig. 13 One 
can see that the spectral weight of the bound state goes 
to zero as it approaches the continuum. 

Next we look at the two-magnon contribution, which 
has only overlap with the even parity operator. In Fig. 114| 
we have plotted different momentum slices of the spectral 
function. At momentum zero, the spectral function is 
identically zero (the ground state is a singlet) and grows 
for small momenta as oc n 2 (cfr. Ref. [ 551 ) . For larger 
momenta, we see that the spectral function gets strongly 
peaked at some value for k. after which the peak again 
disappears. The origin of this resonance is of course the 
formation of the bound state: before it is stable, the 
bound state is already visible in the spectral function as 
a resonance. 

To further confirm this picture, we have plotted the 
maximum of the peak in function of the momentum in 
Fig. [l5j One can see the resonance clearly diverging at 
the point where the bound state reaches stability: from 
that point on the stable bound state contributes a delta 
peak to the spectral function. 

We have also plotted the integrated spectral function in 
Fig. 16 Before the formation of the bound state, we see 
that the sum rules are completely satisfied (up to numer¬ 
ical errors), which shows that the one- and two-particle 
sectors indeed capture the full spectral function, at least 
in this momentum range (see also Ref. l82|) . Again, we 
clearly see the oc n 2 dependence at small momenta. Af¬ 
ter the bound state has formed, however, the two-magnon 
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FIG. 13. The one-particle spectral weights; these appear 
in the spectral functions So/ n (K,uj) as the prefactor of the 
2n5(co — A(ft)) function (where A(ft) is the dispersion rela¬ 
tion of the particle). We have plotted the magnon weights 
w.r.t. to the odd operator (green) and the weight of the 
triplet bound state w.r.t. to the even operator (blue). All 
the other one-particle spectral weights are identically zero. 
These results are in accordance with Ref. l82l Note that the 
one-particle description of the bound state gets worse when 
coming closer to the continuum, so that the calculation of its 
spectral weight loses accuracy in this region. It is nevertheless 
clear that the spectral weight goes to zero as the bound state 
loses stability. 



FIG. 14. The two-particle contribution to the spectral func¬ 
tion So(k,lo) for equally spaced values of the momentum be¬ 
tween k = 0 and k = 7t/2. The ft = 0 curve is not shown as it 
is equal to zero everywhere. Calculations were done at 7 = 2 
with bond dimension D = 32. 


part loses increasing spectral weight to the bound state. 


D. Magnetization process 

Let us now turn on the magnetic field. For SU(2) in¬ 
variant systems, this perturbation does not affect the sin¬ 
glet ground state and induces a Zeeman splitting of the 
elementary magnon excitation. When the magnetic field 
reaches the value of the gap, one of the components of 
the triplet forms a pseudo-condensate (no real conden¬ 


FIG. 15. The maximum of the two-particle contribution to 
the spectral function So (ft, oj) for different momentum slices. 
The full line is a guide to the eye. In the inset we show a close- 
up of the small momentum region, the full line is quadratic 
fit. Calculations were done at 7 = 2 with bond dimension 
D = 32. 



FIG. 16. The integrated spectral function f Alo/2'kS(k,lo) in 
function of the momentum ft (red dots) compared with the 
momentum space correlation function So (ft) (blue line). In 
the inset we plot the (loglO of the) difference between the 
two; values below ICC 2 are not shown. Calculations were 
done at 7 = 2 with bond dimension D = 32. 


sate can form in one dimension); the system undergoes a 
continuous phase transition from a commensurate phase 
with zero magnetization to an incommensurate phase 
with non-zero magnetizatiorpS. 

The physical picture of this condensation can be un¬ 
derstood from the approximate Bethe ansatz that was 
developed in Sec. Ill C Indeed, once it crosses the gap, 
the magnetic field serves as a chemical potential for the 
+1 component of the magnon triplet (the other compo¬ 
nents remain gapped, so we will not consider them in 
our calculations). The information on the magnon dis¬ 
persion relation and the magnon-magnon S matrix we 
have gathered in the previous sections will allow us to 
compute both thermodynamic properties and correlation 
functions for the magnetized chain. 

We start very close to the phase transition, where only 
the momenta around the minimum will be occupied, so 
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magnetic field 



FIG. 17. The magnetization of the ladder (7 = 2) in func¬ 
tion of the applied magnetic field h. The dots are calculated 
with a direct MPS optimization (using an adapted version of 
Ref. [32}, the red line is the free-fermion result [Eq. (411], the 
green one is with the scattering length correction [Eq. @], 
and the blue line is a full approximate Bethe ansatz calcula¬ 
tion. 


that we can approximate them as free fermions. If we 
introduce a characteristic velocity v for the magnon dis¬ 
persion around its minimum as 


A(k) — A + (tO ttmin) ) 

the magnetization (i.e. the density of condensed 
magnons) will be given by 8 ®^ 


m(h) = - h c ). 


7 TV 


(41) 


When more pseudo-momentum levels are filled up, the 
two-particle S matrix will deviate from its limiting value 
of —1 and the free-fermion approximation will no longer 
hold. As a first order correction, we can assume a lin¬ 
ear scattering phase with the scattering length a as the 
slope (and still a quadratic dispersion). From Eq. (32) it 
follows that the correction to the magnetization curve is 
given by 

m(/0 = —(42) 

7 TV 6tt z v z 

a result which was obtained in Ref. [523 by a similar rea¬ 
soning. 

When even higher momenta are occupied these approx¬ 
imations (quadratic dispersion relation, linear scattering 
phase and Galilean invariance) will get worse and only a 
full Bethe ansatz calculation will give the correct magne¬ 
tization curve. In Fig. [IT] we have plotted this. 

Next we look at correlation functions of the magnetized 
ladder. With our methods, we have no direct access to 
these correlation functions, but we can infer their form by 
combining the Luttinger liquid formalism with the ther¬ 
modynamic properties computed from the approximate 
Bethe ansatz. Indeed, since we have seen in Sec. IIV C 


that the S% operator essentially creates a magnon out of 
the vacuum at momentum 7 r and the Sq operator creates 


FIG. 18. The LL parameter in function of the magnetization 
for 7 = 5 (blue), 7 = 2 (red), 7 = 1 (green) and 7 = 1/2 
(magenta). 


a two-magnon state at momentum 0 , we can translate 
the expressions for the Bose gas correlators [Eq. (30)] to 
the magnetized ladder as 




-B,(-l)‘-‘ , C | OS( 2 ! , ”)l: ) 5/ > (43) 


|j _ 1 2K+1/2K 


((S5)AS5)i) = m 2 - 


K 


2tt 2 \ i — i'\ 


A 


cos(27 Tm(i — *')) 


i>\2K 


(44) 


in accordance with Ref. GUI The power-law decay of these 
correlation functions is controlled by the LL parameter 
K. In Fig. [IS] we have plotted K in function of the mag¬ 
netization m for the ladder at different values of 7 . At 
very low magnetization m —> 0 the LL parameter reaches 
the universal value of 1 , but it appears that, beyond this 
limiting value, K(m) changes qualitatively as we vary 7 . 
The same behaviour was observed in Ref. [501 by fitting 
the analytic form of the correlation functions (43) and 
( [44] ) with numerical calculations. 

This behaviour can again be explained by starting with 


the free-fermion limit at very low densities. In Sec. Ill D 
we have shown that the LL parameter equals K = 1 
in this case. The first order correction on this value is 
determined by the magnon-magnon scattering length; in 
first order in m the LL parameters is given bjimi 


K(m) = 1 — 2 am. 


(45) 


In Fig. [l9] we have plotted the scattering length in func¬ 
tion of the interchain coupling 7 . Based on Eq. (45), the 


change of the sign of a confirms the varying qualitative 
behaviour of K{m ) as observed in Fig. 18 and in Ref. [90, 


Finally, we can study the magnetization process at 
finite temperatures using the thermodynamic Bethe 
ansatz. In Fig. 20 we have plotted the magnetization 


curve for different temperatures, showing that the zero- 
temperature square-root dependence around the phase 
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FIG. 19. The scattering length for different values of the 
interchain coupling 7 . 



FIG. 20. The magnetization in function of the magnetic field 
h for three values of the temperature: T = .01A (blue), T = 
.045A (green) and T = .08A (red). 


transition is smoothed out at finite temperature. Note 
that we have included the other components of the 
magnon triplet - they are thermally excited as well in 
a decoupled fashion. In a more correct analysis we would 
have to solve the fully coupled Bethe equations for the 
three components, but this falls outside the scope of this 
paper. 


V. FUTURE DIRECTIONS 


A. Topological excitations and bound states 

In the previous sections we have restricted our frame¬ 
work to the case where we have a unique ground state. 
We can easily extend the framework, however, to situa¬ 
tions where we have symmetry breaking and the elemen¬ 
tary excitations are domain walls rather than localized 
particles. 

Suppose we have a doubly degenerate ground state, 
approximated by two MPS |'F[Ai]) and |SP[Al 2 ])■ The 
obvious ansatz for a domain wall excitation is 


!$.[£]> 


14 


x B 3 


11 -*v 


.m<n 


1 1 


.m>n 


v r |{ s }) > (46) 


i.e. the domain wall interpolates between the two ground 
states. The ansatz has been successfully applied to the 
gapped XXZ model in Ref. 12.11 where the elementary ex¬ 
citations are spinons, and to the Lieb-Liniger model in 
Ref. [55) where topological excitations are elementary. 

Strictly speaking, however, the momentum of the 
ansatz [Eq. (46)] is not well defined: multiplying the ten¬ 
sor A -2 with an arbitrary phase factor A 2 <— A^e 1 ^ shifts 
the momentum with k <— k + (f>. The origin of this ambi¬ 
guity is the fact that one domain wall cannot be properly 
defined when using periodic boundary conditions. 

Physically, however, domain walls should come in 
pairs. The procedure for constructing a scattering state 
of two domain walls is completely analogous as in Sec. [TT] 
For these states the total momentum is well-defined, al¬ 
though the individual momenta can be arbitrarily trans¬ 
ferred between the two domain walls. Scattering states 
of two domain walls are especially relevant as they are 
the first excitations that carry any spectral weight. Con¬ 
sequently, a first non-trivial contribution to dynamical 
correlation functions asks for a solution of the scattering 
problem. 

A second extension of the scattering formalism is to¬ 
wards the study of bound states. As we explained above, 
a bound state should be interpreted as a one-particle ex¬ 
citation and described by a one-particle ansatz. Yet, in 
the case where the bound state becomes very wide - e.g. 
when it is close to a two-particle continuum - the one- 
particle ansatz is not able to capture its delocalized na¬ 
ture. One possible extension consists of working on mul¬ 
tiple MPS tensors at once, leading to the ansatd^U 


In the previous sections we have shown how to varia- 
tionally determine all properties of one- and two-particle 
excitations of generic quantum spin chains. In this last 
section we show how our framework can be extended 
to study domain wall excitations and bound states and 
how to compute spectral functions at finite temperature. 
Since we believe that our work provides a crucial step to¬ 
wards the construction of an effective Fock space of inter¬ 
acting, particle-like excitations, we provide some further 
steps in this direction. Lastly, we reflect shortly on the 
application of our methods to two dimensional systems. 


| *,«[£]> =5> i,en 5>l 

n { s } 


II 1 

_m<n 


x ^S n ,S n + l,...,S n +N 


n v 

_m>N-\-n 


v r IM) ■ 


(47) 


The number of the variational parameters in the big B 
tensor grows exponentially in the number of sites, so that 
we cannot systematically grow the block as the bound 
state gets wider. 
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FIG. 21. Graphical representation of the bound state ansatz. The B tensor of the one-particle ansatz in Fig. [3] is spread over 
more than one site. 


As a more systematic way to study wide bound states, 
we should use the two-particle ansatz ^ to describe 
them. In contrast to a scattering state the energy of 
a bound state is not known from the one-particle disper¬ 
sions, so that we will have to scan a certain energy range 
in search of bound state solutions - of course, with the 
one-particle ansatz we can get a pretty good idea where 
to look. A bound state corresponds to solutions for the 
eigenvalue equation ([8]) with only decaying modes in the 
asymptotic regime. In principle we should even be able to 
find bound state solutions within a continuum of scatter¬ 
ing states (i.e. a stationary bound-state, not a resonance 
within the continuum) by the presence of additional lo¬ 
calized solutions for the scattering problem. 


B. Spectral functions at finite temperature 

At finite temperatures, the thermally excited density 
of excitations already present in the thermal state de¬ 
stroys the perfect coherence of one-particle contributions 
to spectral functions: the delta peaks at zero temperature 
will get smeared out in finite temperature spectral func¬ 
tions. It appears that this thermal broadening depe nds 
heavily on the interactions between the particlej 9 b 92 l , so 
that a full quantum mechanical treatment is needed to 
accurately resolve it. 

At zero temperature the spectral function S(k,ui) can 
be expressed in terms of the spectral weights of the low- 
energy excitations of the system. At finite temperatures, 
this is no longer true as we generally need form factors 
corresponding to states with arbitrarily high energies. 
In gapped integrable systems - where the higher energy 
states can be labelled with a particle number n and have 
an energy of the order nA - the higher-energy form fac¬ 
tors are suppressed with a Boltzmann factor O (e“” A / T ), 
so one can restrict to low-particle form factors at low 
enough temperatures (compared to the gap) 9 ^®2U. 

In this paper we have shown that, even in non- 
integrable systems, the particle picture remains valid at 
low densities (low temperatures), which makes the low- 
temperature expansion in O (e _A//T ) possible for the 
non-integrable case as well (see also Ref. EH for a sim¬ 
ilar expansion for non-integrable systems). So we can 
associate a particle number to higher excitations and we 
can write down the finite temperature expression for the 
spectral function in the Lehmann representation as 

S(k,u>) = E 2t n5 (£({«})-£({/?})-w) 

mn {ct}{[3} 

2 tt 6(K{{a}) - K{{0}) - k) 


B P&({*}) |( m) { am }|0|n, {Pn})\ 2 (48) 


where is a double sum over particle numbers rang¬ 

ing to oo and {a m } is a set of m particle types: the 
states |to, {a m }) can then be identified with the multi- 
particle states in the approximate Bethe ansatz picture 
of Sec. |III C| We can see that, for gapped systems, the 
Boltzmann factor provides a small parameter, so that ex¬ 
citations with many particles only play a limited role at 
low temperatures. In the thermodynamic limit, two diffi¬ 
culties remain: (i) when coming close to the one-particle 
dispersion curve (where the zero-temperature spectral 
function has its S peak divergence) we have to perform 
a resummation in order to take into account an infinite 
number of terms, and (ii) the form factors appearing in 
Eq. (48) can be divergent in the thermodynamic limit. A 


careful analysis shows that both difficulties can be over¬ 
come in the case of integrable (free and interacting) mas¬ 
sive field theorie:P^. Within our framework, it should 
prove possible to calculate finite-temperature spectral 
functions for generic spin chains (non-integrable) and go 
beyond the perturbative approaches of Refs. [Ml and 1571 


C. Effective field theory 


Whereas the approximate Bethe ansatz provides a way 
to construct an effective first-quantized wave function for 
a finite density of excitations, a systematic construction 
of an interacting many-particle model should be formu¬ 
lated in second quantizatiorPSE 27 ^. We introduce mo¬ 
mentum space creation and annihilation operators that 
act on the ground state as 

ct(«)|*[A]> = |*„(«)> 

c a (K) |^[A]) = 0 

and write down an effective interacting theory 


h=yJ ^ a «( k ) c «( k ) c “( k ) 


E 


&K d«i dl\2 
2tT 2-7T 27T 


Va'0' ,a[3 (^b ^1; ^'2) 


x C L( K 1 + k 2 - K)c^,(n)cp(K2)c a (Ki). (49) 


Since we only have explicit access to the operator acting 
on the ground state and not the operator itself, it is a 
priori not clear how to determine the c 7 (k) and c 0 (k) in 
a unique way. Moreover, there seems to be no trivial way 
for imposing the correct commutation relations. Thirdly, 
because these operators will be momentum-dependent, 
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the transition to a local, real-space representation of the 
Fock operators might not be well-defined. The construc¬ 
tion of Wannier states out of the momentum eigenstates 
might provide a good starting point 2 ^, although it is still 
not clear how to find the unique real-space operators 
that are essential for computing the interaction term in 
Eq. (49). 

A different approach can be taken by starting from a 
free theory of particles with generalized statistics that 
match the two-particle S matrix. The following effective 
Hamiltonian 

H o = J a (K)Zl(K)Z a (K) 


indeed captures the low-lying spectrum of the origi¬ 
nal Hamiltonian if the Z a and Z' a are the so-called 
Faddeev-Zamolodchikov (FZ) operators obeying the fol¬ 
lowing commutation relations 

Z a (m)Zp(K2) = S2 S 0 (ki, k 2 )Z s (k2)Z j (ki) 
Z a (Ki) t Z /3 (K 2 ) t = K 2 )Z S (k 2 )^ Z 7 (Ki) t 

Z Q (Ki)Za(K 2 ) t = 2ttS(ki - r 2 )<5 q/ 3 

+ S 5 fi(K 1 ,K 2 )Z S (K 2 yZ J (K 1 ). 

The idea is to look at perturbations of H 0 and express 
them in terms of these FZ operators. Indeed, when ap¬ 
plying a non-commuting perturbation, we could have a 
new Hamiltonian of the form 


H = Hv + Y,i ^ {Mfzi( k)Z p (k) 

a/3 

+ MfZ a (-K)Zp{K) + h.c.) (50) 


where 


M“*(«) = <$ a («)|M|*0(K)> 
M^(K) = mA]\M\r ?a (K,-K)) 


are the particle preserving, resp. particle non-preserving 
parts of the perturbation. For small perturbations, we 
can assume that only small momentum states will be oc¬ 
cupied and that the S matrix is approximately —1. In 
that case, the FZ operators reduce to fermion creation 
and annihilation operators and we can diagonalize the 
Hamiltonian [Eq. (50)] with a Bogoliubov rotation. In 
general, this proves not to be possible® and a more so¬ 
phisticated strategy will have to be developed. 

When studying the time evolution of integrable sys¬ 
tems, the occupation numbers n a (n) = Z^(K)Z a (K) cor¬ 
responding to the FZ operators are integrals of motion^. 
For non-integrable systems this is no longer the case, 
although the observation of so-called prethermalization 
plateaus might point to the fact that they are almost 
preserved. Indeed, the mode occupation numbers n a (n) 
provide a way to distinguish a thejrmal Gibbs ensemble 
from a generalized Gibbs ensemble 100 . Consequently, by 


finding an explicit (real-space) representation of the FZ 
operators we could follow the occupation numbers n a (K) 
through time, also when starting from an interacting the¬ 
ory. 

D. Breaking of integrability and Yang-Baxter 
equation 

Integrable systems possess a number of interrelated 
properties - diffractionless scattering, local conservation 
laws, etc. - that makes them amenable to a number of 
analytical techniques. Once the integrability is broken, 
these techniques are no longer applicable. An important 
question is to what extent the different manifestations of 
integrability survive in an approximate way close to an 
integrable point. 

One simple consistenc y condi tion for integrability is 
the Yang-Baxter equation 1 01 * 102 !, expressing that three- 
particle scattering should be indepedent of the order in 
which it is decomposed into consecutive two-particle pro¬ 
cesses. As such it is a condition on the two-particle S 
matrix. Our methods provide a way to test this condi¬ 
tion for non-integrable systems, and, more specifically, 
to study the breaking of the Ya ng-B axter equation for 
systems close to integrable points 103 . 


E. Higher dimensions 

Matrix product states have a higher-dimensional 
general izat ion called projected entangled-pair states 
(PEPS)P®. Just as in one dimension, it has been shown 
that PEPS are able to capture the ground state proper¬ 
ties of generic two-dimensional quantum spin systemi 103 , 
so it should be able to straightforwardly generalize the 
one-particle ansatz of Eq. ([3| to the PEPS formalism. 
Compared to the MPS setting, however, the computa¬ 
tion of the effective one-particle Hamiltonian is a lot more 
involved, because of the fact that the environment in a 
PEPS contraction is a one-dimensional object itself (com¬ 
pared to the zero-dimensional environment in MPS). 

In Refs. HH and 11061 elementary particle excitations 
in two dimensions were studied by looking at the spec¬ 
trum of the transfer matrix. The next step, i.e. a full 
variational calculation of the effective Hamiltonian ma¬ 
trix, should lead to quantitative estimates of the gap and 
full dis pers ion relations of generic two-dimensional spin 
systems^®. 
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Appendix A: Ground state and one-particle excitations: technical details 


1. Uniform matrix product states 


Consider a one-dimensional quantum spin system with N sites and physical dimension d. In the thermodynamic 
limit (N —> oo) we can define a uniform matrix product state (uMPS) as 


|f[A]> = £ 4 

M=i 


+oo 

n ■ 4 ‘ 

rn =—oo 


^r|{4) , 


(Al) 


which is parametrized by the set of D x D complex matrices A s for s = 1 ,,d or, equivalently, the tensor A £ 
C DxdxD jy can shown that all local expectation values are independent of the .D-dimensional boundary vectors v' L 
and vr if the MPS is injective or pure 31 39 . This is the case if the transfer matrix, which is defined as E = ^2 A s ® A* 
and acts as an operator in a D x D dimensional vector space, has a non-degenerate largest eigenvalue w and the 
corresponding left and right eigenvectors (71 and |r) are full rank when written as semi-positive definite Hermitian 
(D x D) matrices l and r. A proper normalization of the uMPS amounts to rescaling the A tensor as A 3 —> A s so 
the spectral radius of the transfer matrix rescales to unity. Indeed, the norm of the uMPS can be formally computed 
as 


(Ik[A] 14'[A]) = ® v^j E°° g) vrJ oc ( l\r ) = tr (Ir) 


so that a proper rescaling of l and r suffices to fix the norm to unity (the proportionality factor is unimportant as all 
expectation values will contain this same factor). The parametrization of (Al) has a redundancy: the state |\k[A]} 
is invariant under a gauge transformation A s —> G~ 1 A S G with G an invertible matrix. There are different ways for 
fixing this gauge freedom and we will not specify which one to choose. 

The other eigenvalues of the transfer matrix E have significance as well; the second eigenvalue u/ 2 ) determines the 
correlation length £ of the uMPS as 


£ = - 


log(oA 2 )) ’ 


(A2) 


Uniform matrix product states prove to offer a very accurate description of ground states of gapped, translation 
invariant Hamiltonians in the thermodynamic limit. For simplicity’s sake, we will restrict to nearest neighbour 
interactions, so that H = y~L h n ,_ n + 1 . Having found a variationally optimal tensor A for this Hamiltonian (with 
variational energy density eo), we can calculate its variance with respect to the state Ilk [A]) to get an idea of how well 
it approximates the true ground state. This variance scales with the system size, however, so we should define a local 
state error as 


e GS = ^AHgs = (4/[A] | H 2 |tt[A]>, 

where |Z| represents the diverging system size and we have redefined the Hamiltonian as h n>n +i 
simple calculation shows that the local state error is equal to 

ccjs = Tg? EE WA}\ hn,n-\-lhn '|^[-^]) 

1 n n‘ 

= £<*[41* 0 ,lhn' ,n' + 1 \m) 
n' 

+oo 

= 2 x (l\HU £ E n HH\r) + (l\Hf*<£>\r) + |r) + (l\Ji£\r) 


n —0 


where we have used the following notations 


^ A s B t ® C s ' D 1 ' (s't'\ h |st) 


H DE(F)= E A s B t C u ® D“ E" F' x 




(vu'\ h | tu) (s't'\ h |sw) 


J, 


AB 


CD 


= £ A‘S 


•«c*' d‘‘ X 


‘'rijn+l 


- e 0 . A 


^ (vw\ h |st) (s't'\ h \vw) = ^ A s B t <%> C S D* x {s't'\ h 2 |st). 


VW 
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As the transfer matrix has spectral radius 1, the infinite sum does not converge. On every encounter of a geometric 
sum over E, we will seperate it into a disconnected part corresponding to the rank 1 projector Q = |r)(Z| onto its 
eigenspace with eigenvalue 1, and a connected part corresponding to E = E — Q = PE = EP = PEP with P = 1 — Q 
the complementary projector. Since E has a spectral radius smaller than 1, the geometric series over the latter can 
be safely calculated and we obtain 

+oo +oo +oo +oo 

E ^ = E k)(*l + p E E np = ]T |r)(Z| + (1 - E) p 

n—0 n =0 n —0 n —0 

where the extra projector P in the second term is only necessary to ensure the correct results for the n = 0 term, and 
we have introduced the notation 


(1 - E) p = P(1 - E)~ lp = (1 - E)~ lp . (A3) 

(1 —E) p is zero in the eigenspace of (1 —E) with eigenvalue zero, and acts as the inverse of (1 —E) in the complementary 
space. It thus acts as a kind of pseudo-inverse, although we will also use the (■ ■ • ) p notation more generally below as 
(1 — e lK E) p = P{ 1 — e lK ’E)~ 1 P. Now using that (l\Hjfjt\r) = 0 through the redefinition of the Hamiltonian, we can 
conclude that any geometric series of E which has ( l\H^ to its left or H^\r) to its right will have no contribution 
from the disconnected part, and yield a convergent (finite) result. In particular, the result for the “state error density” 
is 


e GS = 2 x - E) p H££\r) + (l\H^/J\r) + (l\H^\r) + (l\J^\r). 


t (A)AA i 


2. The particle ansatz 


The ansatz for an elementary excitation on top of the uMPS ground state, parametrized by the tensor A , is given 
by 


+oo 

*«[£]>= E 

e“E4 

n 

B s " 

n 

n =—oo 

M 

_m<n 


_m>n 


(A4) 


It is the momentum superposition of a localized disturbance, parametrized by the tensor B (same dimensions as A). 
At zero momentum, this excitation lives in the tangent space of the uMPS manifold with fixed bond dimension, at the 
point |’I'[A]) (see Ref. [2S]for more details). The gauge freedom within this manifold has its reflection in the tangent 
plane: the state |$ K [.B]) is invariant under the transformation 

B s -a B s + XA S - e lK A s X 


with X a general (D x D) matrix. The tensors B s = XA S — e lK A s X give rise to so-called null modes. Getting rid 
of them is possible by imposing a gauge fixing condition on the tensors B and introducing a corresponding restricted 
parametrization. Two choices are especially convenient: 

1. Left gauge. We construct the ( qD x D)-matrix A a ,(b,s) = ((A 1 ') s l 1 / 2 )a i t, and find the right null space Vl of L, 
so that LVl = 0. This matrix Vl has dimensions qD x (q — 1 )D and is orthonormalized: Vf Vl = 1. The left 
gauge fixing condition and its reduced parametrization in terms of the ( D(d— 1) x D) matrix X are then given 

by 

(l\E% =0 ->• B l [X] = l-^ 2 V[Xr- 1/2 . 

2. Right gauge. We construct the ( qD x U)-matrix R( a ,s),b = ( rl/,2 (^) s ) Ql b and find the left null space Vr of R, 
so that VrR = 0. This matrix Vr has dimensions (q — 1 )D x qD and is orthonormalized VrV 7 J = 1. The right 
gauge fixing condition and its reduced parametrization in terms of the (D x D(d— 1)) matrix X are then given 

by 

E^\r) =0 ->• B r [X] = r 1 / 2 !^" 1 / 2 . 

The expression for the norm of the state |$ K [S]) is simplified with one of these choices to be just the Euclidian norm 
in terms of the parameters X (up to momentum S factor) 

(^ kI [B l/ r{X')]\^ k [B l/ r{X)]) = 2tt6(k' - k) x tv(X'X) 


(A5) 
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Moreover, with either of these gauge conditions the excitation is orthogonal to the ground state, so that (<f> K [.B] |\I/[A]) = 
2n5(K,)(l\Eg\r) = 0. The overlap of the Hamiltonian between two excited states can be calculated to be (see Ref. 

{<t> K ,[B'}\H\$ K [B]) = 2v5{K-K')[(!\H§A\r) + (l\H£ p |r) + e~ iK {l\H™,\r) + e iK (l\H^ A \r) 

+ (l\E»(l - E) p HH\r ) + - E) p E p ,\r) 

+ e~ iK (l\H£%{l - e~ iK E) p E£,\r) + e~ 2iK (l\H%£(l - e~ iK E) p E£,\r) 

+ e iK {l\H^,(l - e iK E) p E% |r) + e 2iK (l\H^ A (l - e iK E) p E%\r) 

+ e~ iK (l\Hii(l - E) P E P ( 1 - e~ iK E) p E£, |r) 

+ e iK {l\HH{l - E) P E*,( 1 - e tK E) p E%\r) 


(A6) 


for B and B' in the left gauge and 


(* k .[B']\H\* k [B]) 


’/Jl rrR4 


+ e~ iK (l\E%(l - e~ iK E) p H^ A \r) + e~ 2iK (l\E%(l - e~ iK E) p H^,\r) 
+ e iK {l\E^,(l - e iK E) p H AA \r) + e 2iK (l\E^,(l - e iK E) p H^\r) 

+ e~ iK (l\E%(l - e~ iK E) p E^,{ 1 - E) p H$£\r) 

+ e iK {l\E^,(l - e lK E) p E%( 1 - E) p Hi£\r) . 


(A7) 


for a right-gauge fixed B and B'. We have introduced the notation for a “generalized” transfer matrix Eg = 

J2 s a s ®b s . 

Because of the linear parametrization of (A4) in terms of B, variationally optimizing this ansatz can be reformulated 
as an eigenvalue problem 


mm 

x 


(* k [B l/r (X)]\H\$ k [B l/r (X)]) 

(* k [B l/r (X)]\* k [B l/r (X)]) 


Hlp, e ffAl — AN c g jlp A' 


with H lp o ff the Hamiltonian overlap matrix between two excited states (Eqs. (A6) and (A7|) and the effective norm 
matrix Ni Pi eff equal to the identity matrix because of Eq. ( |A5| ). The eigenvalue A is the excitation energy. By repeating 
this procedure for different momenta, we can trace out the excitation spectrum. Note that the interpretation of the 
solutions in terms of one- and multi-particle excitations can be made on the basis of the computation of the variance, 
as explained in Sec. |A4| Indeed, it might very well be that the lowest eigenvalue at a certain momentum corresponds 
to a two-particle scattering state. 


3. One-particle form factors 


The states (A4) provide a variational approximation for the true low-lying excitations of the full Hamiltonian. Their 
overlaps with a local operator acting on the ground state (their spectral weights) provide an important contribution 
to the spectral function 


S(k, oj) = 


+oo 

E 


p+oo 


elf e 


i(ujt— K,n) 


(*o|Ot(t)Oo(0)|tf 0 > 


By inserting a projector on the one-particle subspace, the one-particle contribution can be written as (Ti(k) denotes 
the set of one-particle states at momentum k) 


S(k,u. >)i p = ^2 2TT5(A a (n) - w) ($ K [H a ]| O 0 |^[A]) 


aGri(K) 


The overlap is given by (with B a in the left gauge) 

I Oo MA]) = (l\0£jr) + (l\Oi( 1 - E) p E%Jr) 


where we have again generalized our notation to an “operator transfer matrix” O 


i = E s , t A s mB t < t\d\s). 
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4. Variance of the one-particle ansatz 


If we write the one-particle ansatz as 


* K [B]) = Ye iKn lx(n)> 

with x( n )) = E v l 

11 

B Sn 

11 

n 

« 

_m<n 


_m>n 


where B is in the left gauge, such that the site dependent states are orthonormalized as (x{ n ')\x( n )) — $nn'- The 
variance of the Hamiltonian with respect to this state can be calculated as (we denote A (k) = (<b K [H]| H |<1> K [£?])) 

cex = (M^]| (77 - A(k)) 2 |$ k [H]) 

= E E e ~ iKW (x(n')\H 2 |x(n)> - A(k) 2 ($ k ,[B]\<S> k [B]) 

n n' 

/ +oo 

= 27r<y(«-« / )l E e- w ( X (n')|77 2 |x(0))-A( K ) 2 

\n'=— oo 

Does this expression make sense? First of all, the sum breaks off for high enough n ’, i.e. ( x(n')\H 2 |x(0)} -4 0 
for nl large enough, as we will see later on. The infinite <5-prefactor is also no problem as this is just the norm of 
the momentum superposition state. The contribution (x(0)| H 2 |x(0)) is somewhat more problematic however, as 
it contains an infinite contribution from the ground state error. Therefore, we subtract the (infinite) ground state 
variance from this expression. We get the following 



e E x = 27 tS(k - k') 


+oo 


(x(0)| H 2 - AH GS |x(0)> + Y ( e_W (x( n ’)\H 2 | X (0)> + c.c.) - A{n) 2 


n' — l 


In the calculations it will become clear that this is indeed a finite expression. 

The first contribution. We will first calculate the contribution (x(0)| H 2 — AHgs |x( 0))- We have two inhnite sums 
and one infinite quantity in this contribution, so we have to be precise in our summations. We have 

(X(0)| H 2 - AH gs |x(0)> = Y ((X(0)| h ntn+1 h n \n'+i lx(0)> - (^[A]| h nin+1 h n 'y +1 |v&[A]) (x(0)|x(0))) • 

n,n' 


Every term for n can be calculated individually, making sure that the right number of ground state errors ecs is 
subtracted 


(x(0)|?7 2 — AHqs Ix(0)) 


E [(*1 (#aa(1 - E) p HiiEi + H^Ei + JiiEi + H%>$) (Ei)^~ 3 E B \r) 


n =—oo 


+ mii E {Ei)^- W ^ 2 Hii{Ei)\ n '\- 2 E p \r) 

n '= n -\-2 

+ miiiEi)^ {HiE + E A A HEi + E a EE( 1 - E) p H aa ) \r) - e GS (l\E p \r) 

+ (Z| - E) p H aa E b b + H {a I(a) E b + JaaEE + E^Iab 

+ H AA H BA + HHE b (1 - E) p H AA ) \r) - e GS (l\E B \r) 

+ (Z| (h aa { 1 - E) p H ab + H^) + Jab + + H AB { 1 - E) p H AA ) |r) - CGS (Z|£?f |r) 

+ (l\ (H A i{ 1 - E) p HEi + H^a) + Jba + + h Ea{ 1 - E) p H a £) |r) - e G s(Z|J5f|r) 

I / /1 ( ttAA( i T?\P TtB TjAA i jjAB ttAA , tt{B)AA . jAA 

+ W - b) &B H AA + H AB H AA + H BA(A) + h B J AA 

+ EEHffl£> + EEHiii 1 - EfHii) |r) - e GS (Z|7?I |r) 

+oo 

+ E {v\ - E) p EEE a + HiEEi + HEi ) ( Ei) n ~ 2 Hii\r) 

n —2 

n —2 

+ (l\E§ Y K) n '- 1 <(^) n -"'- 2 <|r) 
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(i\Ei(Eir~ 2 ( 


H 


(A)AA 
AA(A) ' 


T7<A jAA 
P J A J AA 


+ E a H AA{A £ + E£h££{ 1 - E) p H££) |r) - e GS (l\E B \r) 


aa\ 


The infinite sums on the first and last three lines need to be investigated further. Separating all powers of E A into 
connected and disconnected parts, the connected parts will yield finite results. This also enables to interchange the 
sums (with appropriate redefinition of the summation boundaries) in the double sum on the second and second to 
last line, so as to obtain for e.g. the latter 

(l\E B (l-E) p H AA (l-E) p H AA \r). 


The disconnected and potentially diverging contributions that survive in e.g. the last three lines are given by 

+oo n —2 

(l\EE |r) £ [ £ (l\Hii(Ei) n - n '-- 2 Hii\r) + (l\ (h^ + J AA + H AA( ^ + H AA ( 1 - E) p H AA ) \r) - e GS 

n =2 n' — l 

By writing the Yln'=i = Xm'=-oo ~ Sn'=-oo an d substituting n' —n! in the last sum, we obtain 


+oo 

mE\r) J2 [mii( 1 - E) p H AA \r) + (Z| (h { /J AA + J AA + H AA{ ~ A £ + H AA { 1 - E) p H AA ) 

n =2 

- (l\E*\r)(l\H AA (l - E) p ( 1 - E) p H AA \r). 


|r) - £GS 


The terms in the remaining infinite sum exactly cancel thanks to presence of e G s and the finite result of the second 

line is obtained. A similar result is obtained from the disconnected part of the first three lines. Inserting this in the 

complete expression yields 

(x(0)| H 2 - A H GS |x(0)> = (l\ (h aa { 1 - E) P H AA E B + H {A £ AA E b + J££e b + H AA( f£ 

+ OEi + H A A A A E P (1 - E) P H AA ) |r) 

+ (Z| (h aa ( 1 - E) p H ab + Hfif* + JAB + H AB(A) + h ab {1 _ e) p h aa ) |r) 

+ (l\ \h aa { 1 - E) p H ba + H^ BA + J B £ + H BA(A } + H§£( 1 - EfHii) Ir) 

+ W y^AAK 1 ~ tj B ti AA + ti AB ti A A + H BA{A) + ^B J AA 

+ E B H AA ^l + - E) P H AA ) |r) 

- 4 x e GS 

+ (Z| (1 - E) P H AA E A + + h {a £ aa e a 

+Jaa e a + (1 - E) p E p \r) 

+ (i\H AA (i - E) p (h a E + e a a hE£ + e£e§( 1 - e) p h aa ) |r) 

I (j\ rpB (-t -rp\P ( tt(A)AA . rpA jAA , rpA tj-AA(A) 

+ - &) [ H AA(A) + ^A J AA + h A M (A)AA 

+H aa (1 - E) P H AA + e a h aa { 1 - |r) 

+ (Z| (7?^(1 - E) p eEe a + + U|i) (1 - E) p H A i\r) 

- 2 x (Z|iOl - E) p { 1 - E) p H A i\r). 


All other contributions. Next we calculate (y(0)| H 2 |y(l))- No problems with subtracting an infinite amount of 
ground state errors is present here, so we have 

(x(l)l H 2 |x(0)> = (l\ (2 x H A i( 1 - E) p H aa ( 1 - E) p E b a + H {A l AA (\ - E) p E b + j££( 1 - E) p E b 

1 - + 2x^(1- £) P I^f + + Jif) ^ilr) 

+ (l\ (2 x - E) p H ba + + J BA + |r) 

+ (l\ (2 x - E) p E b H aa a + 2 x + H ba(a £ + |r) 

+ 2 X (l\ (H aa ( 1 - £) p £^ + h£%e£ + H ba ) (1 - i;) p ^|r). 
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Analogously, 


<X(2)| H 2 |x(0)} = (Z| (2 x Hii( 1 - E) p HH( 1 - £) P 4 + ff£>£J(l - 

+ J%£( 1 - E) P E B + 1 - ^) p 4 

+ 2 x #4(1 - E) p h£ b + + #^ )AP + 4a) {Ei)Ei\r) 

+ (*| (2 x Hi£( 1 - E) P H B £ + # A *4 + 4a + #4 ( PA ) 4k) 

+ (Z| (2 X #4(i - e) p e b h££ + 2 X h£ b h££ + + h^) 

+ 2 x (i\ (#4( 1 - #) p 4a + h£%e£ + h b £) (#4 + 4(i - #) p #4 


) k) 


and for n > 2 


(X(n)| # 2 |x(0)) = (Z| (2 x #4(1 - #) P #4(1 - #) P 4 + #^(1 - E) p E b 

+ 4i(l - E) p E b + # AA i A j(i - #) p #f 
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+ (z| (2 x #4(i - #) p #4 + # ( ^ } + 41 + #4 ( PA ) 04)"~ 2 ^k) 
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+ 2 X (Z| (#4(1 - ^) P 4a + Hi*E£ + H b £) 
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+ (Eir~ 2 H£i + (E£) n ~*E£( 1 - #) P #4 I |r). 
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tAA 


We can throw everything together in order to obtain 

OO 

^e--"(x(n)|F 2 | X (0)) = 

n=l 

e“4l (2 x #4(1 - #) P #4(1 - #) p + #il )AA (i - #) P 
+4a 0 - ^) P + 4^4(1 - ^) P ) - e-^E)" 1 ^Ir) 

+ e _ift (Z| (2 x #4(1 - #) p #4 + #£J<£> + #££* + 4|) (1 - e-^Ey'Eilr) 

+ e- iK (/| (2 x #4(1 - #) p #4 + H^ ( ab + 4b + |r) 

+ e" 2 -(Z| (2 x #4( 1 - #) p #4 + #£*<$ + 41 + (1 - e-^^-^ilr) 

+ e _iK (Z| (2 x #4(1 - E) P E B H AA + 2 x #44a + H B *£> + k) 

+ e- 2 -(Z| (2 x #4( 1 - e) p h££e b + 2 x h£ b h££ + + #i P)AA ) k) 

+ e- 3lK (Z| (2 x #4( 1 - E) p E b h££ + 2 x #4#4 + H ba £ A 2 + H { a B 22 A ) (1 - e-^r^k) 

+ 2 x e—(Z| 44(1 - E) p E b £ + #44 + #4) (1 - E) p Hii\r) 

+ 2 x e~ 2iK, (l\ 44(1 - #) p #4 + 4a 4 + H%£) (1 - e~ iK E)~ 1 (h££ + 4(1 - E) p h££) |r) 

+ 2 x e- 3i 4| 44(1 - E) P E B £ + h£ b e£ + H B £) (1 - e"* E)~ l H££\r) 

+ 2 x e _4ift (Z| 44(1 - E) P E B £ + h£ b e£ + H B £) (1 - tT™E)~' h££{ 1 - e~ iK E)- 1 E£\r). 

Note that the infinite sum could give rise to one potential divergence coming from the disconnected contribution of 
the last line of (x( n )4 2 l\'(0)) corresponding to 
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However, the first factor is automatically zero if |'I f [^4]) is a variational minimum within the MPS manifold, as it 
corresponds exactly to the directional derivative of the energy expectation value in the direction of B. 


Appendix B: Two-particle excitations: technical details 

In this appendix we give all technical details concerning the two-particle ansatz that was defined as 

+oo M n 

\T(K,uj)) = EE^J 71 ) I Xjj<(n)) 

n =0 j =1 
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n --e- 
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Just as in the case of a one-particle excitation, there is a gauge freedom in this ansatz. We can again choose a left 
or right gauge fixing condition on the tensors Bj , depending on the situation. We will choose to put all B tensors 
in the left gauge fixing condition, which has the consequence that the states | Xj{ n )) are not orthogonal for different 
n (see further). As was argued in the main body, this choice allows for the strictly local term, for which we keep all 
variational parameters, to correct for the inability of the other terms to describe the deformation of the particles as 
they approach. Alternatively, one could choose the left tensor Bj 1 to be in the left gauge and the right tensor Bj 2 in 
the right gauge; this would make the states \Xj( n )) orthogonal for different n. When studying bound states with the 
two-particle ansatz, this might prove to be a better choice. 

Similar to the one-particle case, we can enforce the gauge fixing conditions by implementing an effective parametriza- 
tion in terms of a matrix X with D 2 (d — 1) elements. As we keep all variational freedom in the strictly local term 
|Xj,if(0)), this will correspond to D 2 (d — 1) variational parameters. In the non-local terms |X(ji,j 2 ),if ( n )) we insert 
a basis of left-gauged tensors Bj x and Bj 2 which both describe the (relevant part of the) one-particle spectrum. If 
we have L particles in the system and we need M basis vectors to describe the dispersion of each, we will have 
( L x M) x {L x M) basis states lx(ji,j 2 ),if ( n ))- The gange fixing and normalization conditions on all the B tensors 
can be summarized as 


(1\E B a j = ml 31 = ml 32 - 0 and (l\E%£ \r) = 5 jlJa . 


1. Effective norm matrix 


The effective norm matrix (N e ff)n'j',nj = (Xj'Mn 1 ) |x j,K(n )) has matrix elements 

<XAk(0)IxTk(0)> = 2nS(K - K'){l\E%, \r) = 2nS(K - K')5 hj , 

= 2nS(K - K')(l\E(Eif^E^, |r) 

3 2 

(Xj>,K(n)\Xj,K(n)) = 2 nS(K - K')(l\E(Ei^E^f |r) 

{Xj',K(n')\Xj,K{")) = 2*5{K-K'){l\E%“(Ei) n - 1 E%^Ei) n '- n - 1 E£ \r) 

J 2 


(n r > n). 


2. Effective Hamiltonian matrix 

The effective Hamiltonian matrix (H e f j) n 'j',nj = {Xj' >K (n')\ H \x j,K{n)) has matrix elements 

{Xj',K' (0)| H \ Xj , K m = 2nS(K - K') ' 

m§^ A \r) + |r) + e~ iK (l\H%* |r) + e iK (l\H*»/ A \r) 

+ {1\eZ (1 - E) p H^\r) + - E) P EZ, |r) 
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+ e~ tK {l\H^ (1 - e~ iK E) p e£., |r) + e~ 2iK {l\H% A {l - e~ iK E) p E^, |r) 
+ e iK (1 - e lK E) p E% |r) + e 2iK {l\H^ A (l - e iK E) p E Bj |r) 

+ e~ iK (l\HH(l - E) P E% (1 - e~ iK E) p E£., |r) 

+ e iK (l\HH(l - E) P Ei jt (1 - e iK E) p E% |r)' 

\Xj,k(0)) = 27 r6(K - K')' 


(1\Haa( 1 - E) p E B B [f B i |r) + (l\H2Z E§ |r) + |r) 

3 1 J 2 J 1 J 2 3 2 

+ (l\E% H^ A \r) + (l\Er B i (1 - E) p Hii\r) 

31 3 2 H J 2 

+ e~ iK ((l\Hii(l - E) p E%E^ b „ |r) + (l\H^E^ B \r) + e£ # |r)) 

+ e~ 2lK (l\(H$i( 1 - E) p E B A ’Ei + H^Ei + (1 - e~ iK E) p E^ Bj , \r) 

+ e iK - E) p e£%„ |r) + mil, E B i |r) + (l\H A %_, |r)) 

\ Ji Jo •?! Jo J 2 / 


T^Bj T-iA 


+ e 2iK (1 \(h aa (1 — E) p E b a b + H ab E b + H b ^ b Vl — e iK E) p E Bj |r) 
V 3\ ^2 °1 J 2 J 1 J 2/ J 

(XU'i, 3 ' 2 ),K'{n’)\H 1x^(0)} = 2 t rS(K - K') ' 

(l\(Hii( 1 - E) P E B > , Ei + H AB >, Ei + H B ^ a ) (Ei) n '- 2 E£., |r) 

n — 3 

+ (l\E B i ( £ (E A yH AA (E A ) n '~ 3 ~ i }E B t |r) 

3 1 \ - / ^2 


+ {l\E%, (Eir'~ 2 (Hil, + E a H b ^ a + Ei 

j\ \ 0 2 J 2 


+ e- iK (l\[Hii(l-E) p EZZ +H™*EA +H?g 


P rpBjA jjABj j-iA 


E^,(l-E) p Eii)|r) 

% A ){.Ei) n '~ l Ei |r) 
•?1 / J 2 


/ 71 / jtAA 


Hii(l - E) p E%Ei + H^Ei + H B f) (1 - e~ iK E) p Ei., (Eif E^.Jr) 


+ e iK (l\(Hii{ 1 - £) P ^., + HU , E? + E^ A ) ( E A A y- 2 Ei„ |r) 

\ J 1 / J 2 

+ (Z| - E) P E^, Ei + Hi ^ ( 

n x —1 

^ e ijK (Ei) j - 2 E Bj (Ei) n '- j - 1 Ei. l 


+ e in ' K {Ei) n '- 2 E B 2 + e J ( n,+1)K (Ei) n '- 2 Ei_ ( (1 - e iK E) p E 

3 2 J 2 

(Xo;^),K'(l)| ^ IX(ii A),*:(!)) = - *") " 

(^aa (1 - E) p E B iE B y |r) + E** \r) + mlisi 10 

3 1 J 2 J 1 J 2 J 1 J 2 

+ (l\E B iH B ii\r) + (l\E B iE B y (1 - E) p HH\r) 

3 1 J 2 °1 3 2 


+ e-^ (Z| Hii(l - e)^e;- e;- + H'ii" Ely + e;-;- e^, 


Hii( 1 - E) p E B “E B ” + Hii" E^ + (1 - e~ iK E) p Ei„ E*., |r) 

/ J 1 J 2 


+ e iif (Z| (Eii(l - E) p Ei„ E B B y + HU, E%y + H™%\e%' |r) 

V 3 1 J 2 J X J 2 j 1 J 2 J 


+ e' 2K (l\{Hii{l-E) p Ei,Ei, +Hii„Ei_, + HU B J (1 — e lK Ey E A l E 

\ 3 1 J 2 J X J 2 3 X J 2 J 


(X(j[,j' 2 ),K'(n)\H | XUuj 2 )M n )) = 27T 8(K - K') y 

(i\ (hH( i - E) p E B y (sir- 1 +Hf B y (Eir- 1 + H B ii(Eir~ 2 
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n —3 


+ E2] ^E^E^E^)"" 4 " 3 ^ Ir) 
J 1 / J 2 


*=0 

+ (Z|E^) {E A y~ 2 (h^ + E A H~y2 + E A E^ (1 - EYH aa ) Ir) 

Jl \ J 2 J 2 J 2 / 

(^aa(1 - E) p E^E a _, + H^E A j{ + E^) (E A ) 

rAAn t?\P tt b B t?A , u AB n j^A 


+ e 


E A n E$ \r) 
J 2 


+ (l | (E£*(l - Ey E a 31 E a + H aa ^ 

n—1 

L ^ 2 (^) 


53 (E^ 

— ji 

i=2 


AA 
4V'- 1 


+ e~ inK (Ei) n ~ 2 Ely {Eiy - 1 

H 

+ e- i{n+1)K (Ei) n - 2 El J2 (1 - e~ iK E) p E£ (E^)"- 1 )-^, |r) 

Jl / -^2 

+ e^a|(F^(l-E)^. ; E^ + E^.,Ef" + E^)(e2)"" 2 e£, ^ a |r) 


J B-/ 

J 2 


+ (Z| (E£?(l - E) P E^ ; E a a + E^, E a + H aa 


Jl 
,Bo 


53 e^E^" 2 ^ 1 (E^)""^^, (Ei) 1 '- 1 
i =2 

+ e mK {Ei) n - 2 Ely {Ely- 1 

J 2 

+ e i(n+1)Ar (E^)”- 2 E^., (1 - e iK E) p El jl (E^)" -1 )Ef 3 ' 2 |r) 

J 2 / 

j;,j'),if'(2)| E IXfr^a)) = 2ttS(K - E') ' 

(Z|(e^(1-E) p E^)E^ +HllyE% 2 +Hly^)E A \r) 

+ (Z|E^ (Ef-.f + E B y 2 H AA , A + Ef-E^., (1 - E) P E^) |r) 

Jl \ ■? 2 J 2 J 2 / 


+ e 
+ e 
+ e' 


-iif 


Ji 

rAA I 




T ABi 


#aa( 1 - E) p E^Ely + nyy^Ely + E^Me^E*, |r) 

Jl Jl Jl / J2 

4.4(1 - E) P E^E^ + H AB ^E B y 2 + Ef- 5 -) (1 - 

E^(l - £) P ^., < J1 + , Ef- + H AB y)E n B y |r) 

J 1 Jl Jl / J2 


e~' K E) p E A E A E A „\r) 

Jl J2 


+ (Z| H AA ( 1 - E) P E B E a + E^ E^ + E b 

Jl J 1 : 


y){ 

- Eb2 p Ea 2 + e ”“^ 2 p(l - e iK E)^E P - El 
K'(n + 1)1 H I XUuhhKW) = - *") 


E44 (1 - E) P E|- (E^)- 1 + H^ly {Ely - 1 + e|^(e^- 2 )e^e£ |r) 

Jl Jl Jl / J 2 

n —3 

+ (Zl^l;, 1 (Y.( E A) iH AA{Eiy- i - 3 )E B A i2 E A j , |r) 

Jl i=0 2 

+ (i\Ely (E A y~ 2 (Hfy 2 e a „ + E A H B yyl 

Jl \ J 2 J 2 

+ E A E B y 2 H AA , A + E A E B y 2 E A _, (1 - E) p H aa ) Ir) 

+ e-^(Z| (E^(l - E) p E^E A ]{ + Hl B A ^E A o , + Ef^f) (Ei)"- 2 Ef-E^E4 p |r) 

+ (Z| (E^(l - E) p Ef- E^ + E^ 1 Ei + Ef- A ) ( 


n— 1 

E- 

j=2 


e -i jK ( E Ay- 2 E A t (Eiy-^E^ 2 {Eiy'- n+j ~ 1 E^, l 

Jl J2 
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+ e~ inK (Ei) n - 2 E^ [E A ) n '~ 1 E Bl 

31 3 2 

+ e~ i{n+1)K {Ei) n -‘ 2 E 1 ^ 2 (1 - e~ iK E) p E^„ (Eif^E^., ) |r) 

31 3 2 / 

+ e iK (l\ (Hii( 1 - E) P E^.,Ef« + E^ + H^iEir^E^ |r) 


+ (l | (#ii(l - £) P ££, E A + H%£ E% + H^ a 

V 3 1 Jl Jl 

n'-l 


V e ljK (Eiy- 2 E A jl (Ei) n '- : >- 1 E^_ l (E A ) n ~ n ' +j ~ 1 E A i2 

J=2 

+ e ir! ' K (Ei) n '- 2 Ef J1 (E A ) n ~ 1 E A j2 

32 

+ e i{n '+ 1)K (Ei) n '- 2 E£ , (1 - e iK E) p E A jl (Ei) n ~ l E Bj A\r) 

J 2 / 

*'("')! H IXUurtM 1 )) = ™(K ~ *")" 

(*i fe 4 (i - ^) P ^, 1 ^f 32 +^C 3 , 1 ^f 32 +^f 32 ) (Eir'- 2 E£ | r) 

\ 3 i J! J X / J 2 

+ (l\E^H^ A (Eir'~ 3 E^ \r) 

3 1 J 2 

n'—4 


+ (l\E*“E% 2 J2 (EiTHHiEir'-^E^ |r) 

31 i=0 32 

+ (Z|£7f J > E A 32 (. Eif~ 3 (JT^„ + E A H B y A + E A E B (1 - E) p #ii) |r) 
•?1 \ -?2 -?2 ^2 / 

+ e~ iK (l\(Hii( 1 - E) p E B A ^E B B y + H^E B B y + |r) 

\ Jl •?! 3 1 / 3 2 


+ e 


-22K 


#aa(1 - E) p E b ^E b ’ 2 + HZ’ 1 + <4 


x (1 - e~ iK E) p Et, {Ei) n '- l Ey_, |r) 
•?1 J 2 


+ e‘*(Z| to(l - E) P E B t E B » + Hj&„ E A jl + H B B ;\\E A i2 (Ei) n '~ 3 E£ \ r) 

\ 3 i / J 2 

+ (Z| (Hii(l - E) p Ei, Ei + Hii Ei + 

V 3 1 J X J X 


n —2 


e ijK {Ei) j - 2 E* jl E% 2 {Ei) n '- j - 2 E%, 

3= 2 

+ e i(n '- 1)K {Ei) n '- 3 E A 31 E Bj2 

3 "2 

+ e in ' K (Ei) n '- 2 E Bil E A j2 

32 

+ e i(n'+l)K( E Ajn'-2 E A.' (1 _ e iK E) P E A n E A 12 ) |r) 
H |X(ji ,j 2 ),K{n)) = 2nS(K - K') 


Hii(l - E) p E B yEi + H AB i3 E a + H B yi)(Eir~ 2 E B 32 Ei(Eir'- n - 2 Ei 

3 1 Jl Jl / J 2 


31 i=0 2 

+ (Z|e|J (^)- 2 (tf A f 22 Ei + EiH B % A ) ( Eir'~ n ~ 2 E £ |r) 


n —n —3 


—n —2 ttAA 

h ab., 
J 2 


+ (l\E B y {Eir~ l E B / 2 ( 53 (EiyHii(Eif- n - 3 -*E£_, + (Eif- 

3 1 \ 3 2 

2=0 

+ (Ei)”'-"- 1 Eii A + (Eif'-^Ei., (1 - E) p Eii) |r) 
•?2 °2 / 


+ e _ * K (Z| (E aa (1 — E) p E A 31 E B , +K B “E£., +H B ^ A 

v J 1 -^l ^1 
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x (Eir- 2 E^(Ei) n '~ n El,\r) 


+ (Z| (Hii( 1 - E) P E^Ei + Ei + H^ A ) ( 

n— 1 

( ^ e-* K {Eiy- 2 E£, (Eir-^E^ (Ei)*) (i#)"'""" 1 1 e£„ 

\ z ' 3 \ / 3 2 

J=2 

+ e~ inK {Ei) n - 2 E^ 

3 \ 3 2 

+ e~ i(n+1)K (Ei) n - 2 E f J ' 2 (1 - e-^E^E^., (Ei) n ' ~ x E^ t ) |r) 

■71 -72 / 

+ e^(Z| (1 - E) p Ei E\ 3 » + HU E%' + Hf/X) 

\ ’i -71 / 


x (Eir^E^iEif-^EiJr) 

3 2 


+ (Z| - E) p Ei Ei + HU Ei + H£i A ) ( 

\ 3 1 -7l -7l / \ 

n! —n— 1 

V e ijK (Eiy~ 2 E A Sl (Ei^-'E^ 2 (Eiy'-n-i-^Ei., 

^' -72 

i=2 

+ e i(n ' " n)x (A^)"'- n " 2 J ef J1 ( Ei) n - x Ei 22 

4 

n' — 1 

+ ^ e iiK (Eiy~ 2 E^ 1 (Ei) n '~*~' l Ei (Ei) n ~ n ' +: '~ 1 E 


j=n'— n+1 

+ e in ' K (Ei) n '~ 2 Egi (. Ei) n - l Ei j 


+ e i(n'+l)K( E Ajn'~ 2 E A ' (1 _ A) P A^ n (A^)™" 1 A^ 32 ) |r) 

it / 


3. Asymptotic regime 

The expressions for the effective norm and Hamiltonian matrices above are largely determined by powers of the 
transfer matrices. The power of the transfer matrices behaves as 


(A^)" = \r){l\ + O (e~ n ^) as 


where the correlation length £ of the MPS was defined in Eq. (A2). The asymptotic regime in N e ff and H e g- is reached 
when the corrections can be safely neglected, i.e. n > £ x log(l/e) where e is the allowed error. 

The effective norm matrix reduces to the unit matrix in this regime 

and the effective Hamiltonian matrix is greatly simplified: 

(XUU^K'in )| H \x<ji,h),K(ji)) = 2? tS(K - K') 


S jldi ( 1\{Kb- 1 + H B B ii + Eli (1 - E) P Hii + HH( 1 - E) p El 72 )\r ) 

\ -79 J 2 3 2 3 2 / 


+ 6 jaJ ,(l\(Hii( 1 - E) p Eli + Hili + H B B ii + E B B i (1 - EfHH) |r) 

\ -7^ 3\ 3 1 3 1 / 

OCC iWt),K'{ri + 1)| 77 |X(A,i 2 ),K(n)) = 27rd(A - K') 


5n,jMHii(l-E) p E B A 72 Ei 

\ J 2 


TjABj2 7-iA 

^ AA 


H 


Bj2^ 

AB., 

3 2 


r) 


+ 6 hJ ,e 7K (l | (7?aa(1 - ^ ^ + ^71 17 

\ -7l -7l 3 1 / 

(X(ib4),if'(7)| 7? |X(Ad 2 ),^( n )) = 27T(5(A' - A'') 

<7,4 Gl (Hii(l - E) p E B A 72 Ei + H AB a 72 E a a + Hli A ) (EH n '- n -VEl., \r) 

+ (Z| - A) P A^ Ai + HU Ei + 7^7 ) ( jE ;A)(n'-„-2) jE; B, 1 ^ 
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One can observe that the matrix elements indeed form a repeating row of block matrices, centered around the diagonal 
and exponentially decaying 


(Heff)n'jJj' ,njij 2 — )l H | X(ji,j2),K{n)) 

= (An , -n)j[j^,j 1 j 2 

= O (e“l n as \n — n'\ —> oo. 


4. Two-particle form factors 


Again we start from the spectral function 

+°° /* + oo 


~T /‘-tOO 

S(k, u)= Y / (^ Q \Oi(t)0 0 (0)\^o). 

__ J — OO 


Inserting a projector on the two-particle subspace, the two-particle contribution to this function can be written as 
(r 2 (ft,w) is the set of all two-particle states at that momentum-energy combination) 

S(k,uj) 2p = Y, <T 7 (/e,w)|6o|*o) 

*er 3 (x» 

If we denote the coefficients cP (n) of the two-particle states as 


2T 


«*(«)= 4>caiW + E^"” 
7=1 


such that ^(n) ~ 0 if n > R for some value of R. The overlap appearing in the spectral functions can be calculated 


as 


where 


(M>[A]|6 0 |T(A-, w )) = EE c?(ri) (^[AJI Oo \Xj,K{n)) 

n -0 j 


(vk[A]| do |xj,k(0)) = (l\Of |r) + e* K (l\Oi(l - e iK E) p E% \ r) 

(vk[A]| Oo \ X{h ,h),K) = {l\0 B A p (Eir~'E^ |r) + e iK (l\Oi(l - e iK E) p E f- |r) 

= (l\ (Of 31 + e iK 0%( 1 - e tK E) p E A J1 ) (Ei) n ~ 1 E^ ia |r). 


We have 


(T[A]| Oo |T(A»> = Y ^(0) Ul\0%\r)+e* K (l\Oi(l - e iK E) p E%\r) 


+ E E c wi 2) W(«l (O^ 1 + e lK Oi (1 - e lK E) p E B A -) (E^E^ |r) 

n=lji,j 2 

2r 

+ Y E ^ lda) (l\ (Af- + e iK Oi (1 - e iK E) p E A jl ^j (1 - e iK i n E) p E A i2 |r). 


7 =1 1102 


Appendix C: Proof of equation ( |19| ) 

Let us start with the polynomial eigenvalue equation for the asymptotic solutions of the scattering problem, Eq. (13) 


+M 

Y 

m——M 


V = UJV. 
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For a given /i € C, this equation is an ordinary eigenvalue problem with eigenvalue oj and eigenvector v. Given 
the property A)' n = A_ m , we are only assured of a Hermitian eigenvalue problem (and thus of real eigenvalues w) if 
/.t* = /it , i.e. if /.t is on the unit circle. So for any /x 7 = exp(tK), let there be eigenvalues w 7 (k) and corresponding 
normalized eigenvectors v 1 (n). The functions w 7 (k) and Vj(k) are assumed to be smooth such that at least the first 
derivatives are well defined. By taking the derivative of the eigenvalue equation with respect to k we obtain 


+M 

E 

n=-M 


imA n 


i ( k ) ■ 


+M 

E 

n=-M 


A p*‘ 

r 'm 


^( K ) = ^( K k( K ) +W7 ( K )^(.). 


By multiplying this equation with vy(ny and using the normalization vy (k)1v 7 (k) = <5 7 ' 7 , we tobtain the following 
relation for later use 


-\-M , 

E imvy{n)^ A m e IKm v 7 (/v) = (5 7 ' 7 -^(/v). (Cl) 

m=-M K 

Now consider a two particle eigenstate |T(A', w)), which has the asymptotic form 

2r 

c(K,u) = y^g 7 e^ n v 7 . 

7=1 

We can introduce the projectors (we will omit all dependencies on the total momentum K) 

R L n L„ 

p R = E E IX 7 '( n )> (Xj(n) I and = E E M n )> (Xj(n)\ 

n—0 j—1 n>Rj =1 


so that we have 


<T(c)| P^HP R |T(a,)} = (T(w)| PrHP^ |T(w)> 

upon the condition that |Y(w)) is an eigenstate. If we choose R > M, we can insert the asymptotic form for the 
effective Hamiltonian 


R 


E yy c(n , ) t A„_ n 'c(?r) - c(n) t A n /_„c(?r') = 0. 


n—0 n'>R 

Since A m = 0 for \m\ > M, this allows to restrict the summations and rewrite this equality as 

M R 


T, T, c(n) f A m c(n + 


m) 


.m =1 n=i?+l —m 


= 0 . 


We can insert the asymptotic form for c(n) to obtain for the “diagonal terms” where 7 has the same value for both 
sums 


' 2r m 

Ei^i 2 E 

.7=1 m—1 n—R-\-l—m 


R 

E v\k mV y^ m 


■ 2 r m 

E l^l 2 E ™-v\A m v~ f e ZK 'i' 

.7=1 771—1 

2T M 

= - E y 2 E (* TO « 7 A m« 7 e lKTm - 

7=1 m=l 

2 r M 

= -y\ q y y imv\A m v 1 e l ^ m 

7=1 m=—M 

2T 

= -Ew a ^<%> 

7=1 


im-utA m t i; 7 e-^ ro 


and this expression has to equal zero if we can show that the contribution of all “non-diagonal terms” (7 y 7 in 
the two sums) vanish. We look at a single contribution with 7 / 7 ' and the corresponding term with 7 and 7 ' 
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interchanged, and first assume ^ 
m ft 


M ft 

<lY<hY. J2 v\,A m v^ m e l ^-^ n + q ; qy ]T 51 ^A m1 ; r e ZK v m e l ^'- K ^ n 

m—1 n=i?+l —m m= 1 77,=ft+l —m 

e i(K y -Ky)(R+ 1) JH. 


E 4 A ™<7 


i(rey-re 7 )(fl+l) J£ 

e ,. , £*? E '4 A ™«V - e“’”) 


ra=l 


l 


2 sin((re 7 — ny)/2) 


3? 


M 


3 i(« 7 -Ky)(fl+i/ 2) 9 * /?7 ^5 v t,A m t ; 7 (e ire ' 1 ' m - e iK -y /m ) 

m= 1 

M 

- e l(K v- K -r)(^+i/2) g * g7 , ^ „t A m « 7 , ( e “-r ,m - e^ m ) 


2 sin((K 7 — ny)/2) 


3? 


M 

e i(« T -K 7 ,)(fl+l/ 2) 9 * /(?7 £ V t /AmV 


m=l 


7 \ 

M 


2 sin ((/« 7 — Ky)/ 2 ) 


3 ? 


- e i(K 7-v)( iJ + 1 / 2 ) g ;,g 7 55 ( e “ iKym - 

m—1 
M 

3 i(« 7 -K 7 ,)(fl+ i/2) 9 *,g 7 55 A m w 7 (e iK i' m - e if V m ) 




n.=-M 


Note that we are missing the term for m = 0, but that this term is zero anyway because of the factor (e lK -i m — 

Finally noting that J2m=-M A m v^ m = w « 7 and £ m= _ M VyA m e lf V rra = u«y, it is clear that both contributions 
cancel and the total expression evaluates to zero. Finally, we consider the case that « 7 = Ky = n. We obtain 


MR. MR 

+ 9 7 55 55 v! f ,A m v 7 e lKm + q * qy ^2 55 v^A m vye %Km 


m=l n=ft+l—m 


771—1 77,=ft+1 — 777, 

M 


M 


+ 97 55 ™^ 7 'A m v 7 e* Km + g*gy 55 mv* A m t> 7 /e®‘ 


771=1 

M 


771=1 


9 7 '9 7 


55 mv^Am^e 1, 


m— — M 


In the last line, we replaced the second term of the line before by the negative of its complex conjugate, since we 
are taking the imaginary part of the whole expression anyway. Using that and vy correspond to some v 7 (k) and 
vy(n) with different 7 7 ^ j' but equal uj 7 (k) = cjy(K), we can employ Eq. (Cl) to conclude that this term is zero. 


Appendix D: Mpller operators, the S matrix and scattering states in one dimension 

In this appendix we will translate some basic notions of single particle scattering theory from an external potential 52 
to the one-dimensional case where we have different types of particles with general dispersion relations. The two- 
particle scattering in the many body Hilbert space considered in this manuscript can be mapped to this setting by 
taking out the conservation of total momentum and only looking at the relative wave function, which is encoded in 
the coefficients o’(n). For the remainder of this section, we assume to have a Hilbert space spanned by states {+, j)} 
where a: is a spatial coordinate that can be discrete (x £ Z) or continuous (x £ R) and j = 1 ,L labels different 
internal levels at every position (corresponding to different particle types). We assume we have some Hamiltonian H, 
which can be written as the sum of a free part Hq and a potential V. The free Hamiltonian is translation invariant 
((x',j'\Ho\x,j) = (A x -x')j',j w ith Aj, = (A_ x )t) and also assumed to be short-ranged (A x _ x / = 0 for \x — x'\ > M). 
The potential is centered around x = 0 and goes to zero quickly, e.g. (x\f\V\x,j) = 0 for |x| > M + N or 
\x'\ > M + N. The free Hamiltionian is diagonalized in momentum space and describes the free propagation of 
a number of types of particles a = 1, ...IV with eigenvalues (dispersion relations) E a (p). Indeed, by using the 
momentum states | p,j) = f dxe lpx \x,j) (an integral over x should be read as a sum for the discrete case), the free 
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Hamiltonian Hq is brought into block-diagonal form: 

ip,j\H 0 \p',j') = 2tt S(p - p'){A{p)) jtjl ( D1 ) 

where the Lx L Hermitian matrix A (p) = f dxe lpx A r is an analytic function of p (since A x vanishes for \x\ > M + N). 
Its eigenvalues E a ( p ) and corresponding eigenvectors v a ( p ) define the spectrum of Ho . Also note for further reference 
the relation 


Mp) f ~^(p)v a (p) = (p)6 a ,p. (D2) 

We will henceforth denote the eigenvalues of the free Hamiltonian Ho as E(p a ) and the corresponding eigenstates as 
| p a ) with coordinate representation (x,j\p a ) = v 3 a (p)e lpx . By choosing v a {py vp(p) = 6 a ,p, the eigenstates \p a ) of Ho 
are normalized as {p'p\p a ) = 27 x5(jt — p)5 a jj and span the whole Hilbert space 

1 = j ^ Ip«> (p«i • 

The range of p determines whether we are dealing with a discrete or continuous system, and will not be specified. In 
order to describe scattering experiments, one should build wave packets from these momentum eigenstates 

10a) = J ^ HP) I Pa) ■ 

Typically, we will be interested in wave packets 0(p) that are strongly centered around some momentum po , so that 
it makes sense to express scattering amplitudes (S matrix elements) in the basis of momentum eigenstates. 

Let U{t) and Uo{t) denote the unitary evolution associated to respectively H and Hq- We now want to describe 
some orbit U{t) | ip), which has an in-aymptote and an out-asymptote in the following sense 

U(t) \ip) —> U°(t) IV'in) as t —> — oo 

U(t) I'll)) ->• U°(t) IV’out) as t ^ Too. 

For given jr/'in) or l^out), one can try to define 

10) = lim U(tfU°(t) |0 in ) = fl+ |0in) 

t — y — oo 

10) = lim C/(t) t C/°(t) |0 ou t) = |0out) 

t —>+OO 

with 0± the Mpller operators. The existence of these limits, and thus of the Moller operators, can be proven by 
studying wave packets and linear combinations thereof. For a quadratic dispersion relation, the dispersive behavior of 
the wave packet is often sufficient to guarantee convergence. Since we are studying general dispersion relations E a (p), 
a sufficient condition can be obtained by restricting to wave packets centered around momenta po with non-zero group 
velocity dE a /dp ^ 0. As the limit of unitary operators, the Moller operators fi± are isometries. Finally, we need the 
condition of asymptotic completeness (which is often harder to prove) to ensure that the range of and is the 
same: they map every state to the space of scattering states and satisfy the intertwining relations 

HQ± = n±H 0 . 

The scattering operator or S matrix can then be defined as the operator mapping the in-asymptote to the out- 
asymptote 


|0out) = |0i n ) = S |0in) S = fit £l+. 

One can easily show that the free Hamiltonian commutes with S so it makes sense to represent the S matrix in the 
basis of free momentum states 


(Qis\S I Pa) = 2Tr5(E(qp) - E(p a )) x S qfj<Pi 


If asymptotic completeness is obeyed the S matrix is unitary, which can be expressed in the momentum basis as 

(< 7 / 3 1 S^S\p a ) = 2ir5(p a - qp)5 a 0. (D3) 
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We can translate this condition to the matrix elements S qf)tPa as 
(9/3 I S^S I P a ) = Y^ j Y { q p\S ] M (*71 S I Vex) 

= E / Y^^ E(Pa) - ~ E ^))Sr^S r ^ a 


= E 


dr — 

C C 

2 ^ ° r n ( ip 0r iiPc l 


E 

vp^^Mpo.) 


dE 

dp 


{Pa 


-1 


27T<5(p a / - r 7 ) 2nS(E(qp ) - U(r 7 )) 


E 5 

kr-,eA(p a ) 


r-nqpSr-nPa. 


dE 

dp 


( 3 - 7 ) 


= (S+S) 


9/3Pa 


x 2nS{E(qp) - E(p a )) 
1/2 



dE 

1/2 

dE 

X 


2ir5(E(qp) - E(p a )) 

“T Pa 
dp 


where ri(p a ) is the set of momenta {qp} such that E(q f 3 ) = E(p a ), and we have defined the matrix elements of S as 

- 1/2 

(D4) 



d E, x 

- 1/2 

d£, , 

0 Q/3,P a — 

*<«■> 

^9/3 

-p(Pa) 

dp 


Unitariness of the S matrix, Eq. (D3), implies that S q;i should be a unitary matrix. 

There are different ways to calculate these S matrix elements; one way is to construct the stationary scattering 
states, i.e. the eigenstates of the full Hamiltonian H = H 0 + V. One first introduces the Green’s operators as 

G°(z) = ( z-H °)- 1 
G{z) = {z~H)-\ 

which are related through the relation 

G{z) =G°(z) + G°{z)VG(z) 

= G°(z) + G{z)VG°(z). 

The T operator is defined as 

T(z) = V + VG{z)V 

for which we can easily derive the Lippman-Schwinger equation 110 

T{z) = V + VG°{z)T{z), 

and the equations 

G°(z)T(z) = (G°(z) + G°(z)VG(z))V = G(z)V 
T(z)G°(z ) = V(G°(z) + G(z)VG°(z)) = VG(z). 

The Lippman-Schwinger equation can be rewritten as an integral equation for the matrix elements of T(z) 

M T(z) | p a ) = (qp\V\p a ) + Y^ J ^ {M ^^T(z) | p a ). 

One can derive a related equation for the Mpller operators 

U+|</>)= lim U(t)^U°(t)\4) 

t—t — OO 


r° 

= \(f>) — i / drt/(r)^Uf7°(r) \<f>) 

J —OO 

dTe CT U{T)WU°(T) | 


' —OO 
r 0 


= 0) -l 


= l^)-iE/ Y / dre£T?7 ( r ) tFt/0 ( r ) I PvUPo 


E 


dp 
27r 


G{E(jp a ) + id) V | p a ) ( p 0 
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where we have introduced the time-dependent damping factor to the potential V —> Ve £t , which is allowed for t —>• 0 
according to the adiabatic theorem. The S matrix 

(9/3 1 s I Pa) = {qp\ I p a ) = lim (qp\ | p a ) 

t—>oo 

can be worked out by writing it as the integral of its derivative 

/•OO 

( 9 / 31 5 | p a ) = (q p \p a ) -i dt{q p \ (e iHot V e ~ 2iHt e iH o t + e iH 0 t e - 2 im Vei H 0t) ^ 

Jo 

= (9/3K) -ilimj™dt(qp\ (v e ^ E ^'> +E ^ )+ie - 2H ^ t + e i(E(^)+E( Pa )+ie-2H)ty^ ^ 

= (qp\p a ) + \ 1™ ( 9 / 3 | (vG ^(E(p a ) + E(qp)) + iej + G (^(E(j> a ) + E{q p )) + iej V^J \p a ) 

= (9/3 b a) + lim 


(qp\T[ ~(E(p a ) + E(qp)) + it ) | p a ) 


e ->0 \E(qp) - E(p a ) + it E{p a ) - E(qp) + it 
= 2ir6(qp - p a )5p a - 2ir6(E(qp) - E{p a )) i (qp\T(E(p a ) + iO) \ p a ). 

The off-diagonal elements of the S matrix are given by the on-shell T-matrix elements. We define the amplitudes / 

- 1/2 



d E 

- 1/2 

d E 

/(9/3 t-Pa) = -i 


(qp\T(E(p a ) + iO) 1 p a ) 



which are the off-diagonal elements of S as defined in the unitary matrix (D4). 

We can now define the scattering states 

\p a ±) = \Pa) , H | p a ±) = E{p a ) | p a ±) 

which, through the Lippmann-Schwinger equation for the Mpller operators, obey the relation 
\p a ±) = I p a ) + G(E(p a ) ± iO)V \p a ) = \p a ) + G°{E(p a ) ± iO)V \p a ±) ■ 

Another important relation is 

(qp\T(E{p a ) ± iO) | p a ) = ( 9/3 1 (V + VG(E(p a ) ± *0)V) \p a ) = (qp\V \p a ±). 

An explicit expression for the asymptotic wave functions of the scattering states can thus be obtained: 

(x,j\Pa+) = ( X,j\p a } +Y1 j dx ' ( X ’j\G°(E{p a ) + iO)\x',j') (x',j'\V\p a +) 


= e ,p “V>)+ f dx'^2{x,j\— —— '- —\x',j') {x',j'\V\p a +). 

J “T E(p a ) - H 0 + i 0 


(D5) 


(D6) 


Since we know the exact eigenvalues and eigenvectors of Hq , we will now first introduce a resolution of the identity 


which brings the Green’s function in block diagonal form 
{x,j\p a +) =e lpx v J a (p) + /dz^/gf 


271 \E(p a ) - A (q) + iO 


Jq(x-x') I I ■/ 


j \V\p a +) 


with the matrix A (q) an analytic function of q , as dehned at the beginning of this section. The integral over q can 
be calculated with the residue theorem. For continuous systems, where q ranges over the real axis, we will have to 
close the contour in the upper or lower half plane depending on the whether x — x' > 0 or x — x' < 0. A first set of 
poles will be close to the real axis and c an be obtained from the eigenvalue decomposition of A (q). Together with the 
analytic dependence on q and Eq. (D2), we obtain 


d E. 


A (q ± *0) = ( E (<lp) ± *0 —( 9/3 ) ) v 0 {q)v p {qy. 


(D7) 
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We should therefore separate the set A(p a ) of all solutions qp for which E(qp) = E(p a ) into two parts A ± (p a ) 
corresponding to solutions for which the energy derivative p^(qp) is positive (+) or negative (—). We then find a 

first set of poles of )_ 1 A ( g )+io ) which are of the form qp + iO for qp £ A + (p a ) and of the form qp — *0 for 

qp £ A~(jp a ). The corresponding residues are given by 

Um {q~ {qp±i0)) 

q->qp±i 0 
qfieA ± (p c ,) 

Aside from those solutions, there could be other solutions q = *A 7 further away from the real axis (SKA ^ 0). 
These correspond to values of A where the analytically continued (but non-hermitian) matrix A(zA) has a real eigen¬ 
value £ 7 (zA) = E(iX-y) that equals E(p a ); we denote the corresponding left and right eigenvectors as u> 7 (A) t and 
Wry(X) (wich will in general not be related by hermitian conjugation). The corresponding residue is then given by 

— ^{X 1 )wi 1 w 1 e~ x ^ ( ' x ~ x " > or more generally — (A 7 )Pj,j'(*A 7 )e _A ' 1 '^ x_a: ) with P(zA 7 ) the corresponding eigenspace 

projector. 

Let us now return to the evaluation of the integral over q. Depending on the sign of x — x', we will close the contour 
in the upper or lower half plane and pick up the contributions of the poles in those respective domains. Since we also 
have an integral over x', it seems we will need to split this into the two regions x < x' and x > x'. However, we can 
make use of the locality of the potential to conclude that (x',j'\V\p a +) is only nonzero for \x’\ < M + N. Thus, if 
|cc| > M + N, then x — x’ will have a fixed sign throughout the integral over x’. For e.g. x — x’ , we will need to sum 
up the contributions of all the poles in the upper half plane, corresponding to qp + iO for pp £ A + (p a ) and all iX 1 
with 3?A 7 > 0. The latter contributions will actually vanish if we now take the limit x —> oo. We can then write the 
asymptotic wave function as 


1 


E(p a ) - A (q) + iO 


iq(x-x') _ 


3,3 


- ( f <*> 


v°p(p)v j p{p)e^ x - x '\ 


{x,j\p a +) « v J a {p)e iPaX - i 


S wG A-(p Q )^(9)e^ a 


Y, q ^A+( Pa ) v p(qy qfiX f(??) 


d E 
dp 

d E 
dp 


(qp) 


(qp\V\p±) 

(qp\V\p±) 


X 

X 


and with Eq. (D5) 


(x,j\Pa+) « v 3 a {p)e WaX - i 




En 


\(q)e iq ^ 




(qp\T(E(p a ) + iO) | p a ) 


-l 


—00 
+oo 


(qp\T(E(p a )+ iO)\p a ) x 


X -A — OO 

Too 


^qp^A+{p a ) U P 

The coefficients that appear are the amplitudes that were defined earlier, so we have the nice final result 

- 1/2 


(x,j\p a +) = v 3 a {p)e Wa 


E 

qpeA ± ( Pa ) 


d E, s 

1/2 

d E, , 


f(qp ±- Pa) 

*>> 


P*9/3 X 3 


v°p(q), 


±00. 


(D8) 


For discrete systems, we can proceed in a similar way. Momentum integrals now range from 0 to 2tt and we auto¬ 
matically obtain a contour around the unit circle in the complex plane by going to the complex variable /i = e lq (for 
x — x' > 0) or p = e~ lq (for x — x' < 0). The derivation then follows analogously. 













